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region near the focus, thus indicating the limits of applicability of a small 
perturbation solution. An analysis of the forbidden region should include 
radiative heat transfer effects. 
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INTRODUCTION 


1. Literature Survey 

The prospect of nuclear fusion recently has stimulated a great deal 
of interest in the production of clean, very high temperature plasmas. 
Focusing the beam of a Q-switched laser in a gaseous medium creates a 
microexplosion and an extremely hot plasma is formed. A typical laser 
pulse has a peak power of 100 MW and a width at half power of 20 ns. The 
necessary energy concentration for breakdown corresponds to electric fields 
of the order 10^-10^ V/cm; then gases normally transparent to laser light 
become ionized and hence absorbing. 

The first reported such experiments are by Meyerand and Haught [1,2] 
and Damon and Tomlinson [3] in 1963. Since then a great many investigators 
have approached the problem from both experimental and theoretical points 
of view. The whole phenomenon is usually divided into three distinct 
phases: 1) breakdown mechanism 2) expansion of plasma during the laser 

pulse 3) decay of plasma after the pulse. 

Experimental evidence is considerable and only a very brief survey 
will be attempted, more extensive descriptions can be found in the litera- 
ture and a reference list is given in a very complete bibliographical 
review by C. DeMichelis [4] in 1969. Experiments concerning the first 
phase measure the threshold energy density as a function of various 
physical parameters: gas pressure [5,6,7], frequency of laser radiation 

[8,9], focal length of lens and hence energy losses by diffusion from the 
focal volume [6,7] and even space and time variations of the laser beam. 

The results show the threshold to decrease with increasing pressure with 
a minimum in the very high pressure range (5,000 psi), to exhibit a peak 
as function of frequency, to increase with decreasing focal length and 
to be fairly insensitive to variations of the laser beam. A large 
selection of gases were used: air, rare gases and some metal vapors. 

The most striking feature of the second phase is the speed at which 
the front moves against the laser beam, typically 10-10 cm/s. This 



motion has been recorded using techniques such as streak photograph (first 
by Ramsden and Davies [10] then various investigators [11,12,13]), Doppler 
shift of the radiation scattered by the plasma [10,14,15,16], Schlieren 
photographs [17] and scanning photographs [14,15,16]. 

All the experiments report longitudinal and lateral velocities as 
function of time. A note must be made about the experiments of Veyrie and 
Floux in 1968 [18] and Korobkin and al. in 1968 also [16] which show dis- 
continuous motion of the front in direction of the lens suggesting a 
possibility of multiple breakdown. This effect has been explained by Evans 
and Grey Morgan [19,20]: primary spherical aberration of simple lenses 
causes peaks of intensity along the beam axis thus creating several 
possible regions of breakdown; the spacing between these regions is shown 
to agree very well with Korobkin 's experiment. In fact Veyrie and Floux 
[18] reported that multiple breakdown did not occur when using a lens 
perfectly corrected for spherical aberrations. Other physical quantities 
pertaining to the second phase have been measured: electron density by 
interferometry [21,22] and temperature from X-ray emission of the plasma 
[15,23,24,13]. 

Experiments have been also performed during the third phase. Observed 
were the plasma motion [10,12,25], electron density from microwave trans- 
mission and reflection in a report by Askarayan and al. [26] (they also 
point out the existence of an ionized region ahead of the front) , and 
temperature from spectroscopic studies [27,28,14]. 

The three phases have been subject to theoretical studies. The break- 
down itself is generally described in terms of multiphoton ionization and 
cascade ionization. Since the probability of a free electron to occur 
naturally in the time interval considered (20 ns) is very low, the first 
electrons are believed to be created by a multiphoton ionization process in 
which several photons are simultaneously absorbed, their combined energy 
being sufficient to knock off an electron from a neutral atom. The 
probability of such an occurrence has been calculated and hence the 
threshold fields [29,30,31]. Except at very low pressure, order of 
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magnitude agreement with experiments is obtained only if the computation is 
limited to the creation of a single electron in the focal volume thereby 
indicating that another process takes over to achieve the plasma electron 
density (typically 10 cm ), namely cascade ionization. Free electrons 
absorb energy from photons by inverse bremsstrahlung until they are 
sufficiently energetic to knock off a new electron by collision with a 
neutral atom giving two free electrons [32,33,34]. A cascade develops 
according to N = N o e c ^ T with t: cascade time constant depending on the 
electric field and breakdown occurs if N reaches a critical value during 
the laser pulse. It should be noted that neither of these two processes 
can account for the experimentally observed frequency dependence of the 
threshold field. 

Before turning with somewhat more detail to the theoretical investi- 
gations of the second phase, it may be noted that the last phase has been 
described in terms of blast wave theory [35]. Panarella and Savic [36] 
developed a perturbation theory from a spherical blast wave assuming 
locally radial flow: the shape initially oblong due to asymmetric energy 
addition evolves to a spherical blast wave. 

The second phase, expansion of the plasma under the influence of the 
laser pulse has been theoretically described in terms of three quasi one- 
dimensional mechanisms: a radiation supported detonation wave, a breakdown 
wave and a radiation transport wave. 

The radiation supported detonation mechanism, first proposed by Ramsden 
and Savic in 1964 [35] assumes the wave front to be a shock wave followed 
by a layer in which the laser radiation is absorbed. Thermal reradiation 
from the hot plasma is neglected. Using quasi one-dimensional conserva- 
tion laws across both the shock and the absorbing layer and postulating 
Chapman-Jouguet (normal detonation) conditions behind the discontinuity 
Ramsden and Savic arrived at 

D = [2 (Y 2 -l) ^ ] 1/3 (1-1) 

p o 

where D is the wave velocity, J Q the heat flux and P 0 the density of the 
gas ahead of the wave. Taking into account the conical character of the 
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2 

focused laser beam, J Q ~l/r and Integrating (I— 1) with respect to time 

3/5 

they showed that the radius of the detonation front grows as t for a 
constant power laser and hence 

D = K (V 3 f 2 '5 < M > 

e c 

where K is a constant and P L the laser power. Raizer in 1965 [37] extended 
the model to general hydrodynamic discontinuities (not necessarily Chapman- 
Jouguet) and attempted to include the effect of lateral expansion of the 
plasma. In a not very clear paper by Champetier in 1965 [38] the purely 
one-dimensional plane motion of the gas behind the detonation wave is 
investigated: an expansion wave follows the front whereas an entropy line 
and a shock wave travel in the direction opposite to the lens. The effect 
of time dependence of the laser pulse was considered by Daiber and Thompson 
in 1967 [12]: in the case of a Gaussian pulse shape the time exponent of 
the radius of the detonation is shown to be -3/5 when initial breakdown 
takes place near the peak of the Gaussian and larger for breakdown 
occurring at earlier instants. They also developed a model for gases 
(like hydrogen and deuterium at pressure less than 3 atm) almost trans- 
parent to laser radiation. A last improvement was made by Key in 1969 
[39]: ionization was included in the gas law, the net effect being to lower 
the wave velocity and gas temperature for the same heat flux. 

Champetier and al. in 1968 [40] and Wilson and Turcotte in 1970 [41a] 
independently studied the flow behind a spherical laser-driven detonation 
for constant power addition. Both used a self similar analysis of the kind 
Sedov [42] first proposed for regular blast wave; the hydrodynamic equations 
were integrated in a temperature-velocity plane starting from the saddle 
point corresponding to the origin. It was found that the wave front is not 
Chapman-Jouguet but rather overdriven (i.e. the gas behind the detonation 
is subsonic with respect to the front). Wilson [41b] also considered the 
case of linearly increasing power as well as plane and cylindrical 
geometries. 

The breakdown wave mechanism is based on the idea that during the 
increasing part of the laser pulse breakdown conditions which were satified 
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at time t Q , at the focus assumed to be of cross sectional area A Q> will 
be met further up the beam at location r of larger cross section A at a 
later time t. In 1965 Raizer [37] postulating that breakdown occurs by 
cascade ionization from an original number of free electrons present in 
the cold gas proposed as breakdown criterion that the electron density 
reaches a certain critical value. If the cascade time constant t is 
simply inversely proportional to the heat flux this criterion yields a 
breakdown wave motion characterized by the time integral from zero, 
beginning of the pulse to t of the heat flux at location r being a 
constant. For the rising part of a triangular pulse of maximum power P max 
at time the front velocity is obtained 


D = K (15 ax ) — 

tmax tn ° 


(1-3) 


where K is a constant function only of the electron cascade and the critical 
electron density; a is half the divergence angle of the laser beam. Alcock 
and al. in 1968 [43] remarked that the free electrons necessary to initiate 
the cascade could be produced by precursor ionization due to the thermal 
radiation emitted by the hot plasma. They modified Raizer's model by 
saying that the initial electrons exist only at time t^ just ahead of the 
wave; the breakdown motion is then described by the time integral from t^ 
to t of the heat flux being constant. Alcock and al. [43] further assumed 
(without any theoretical basis) that the time lapse t - t^ between pre- 
ionization and actual breakdown at any station is a constant t^. Under 
these hypotheses the wave velocity foi^a tr ^ an 8 u ^- ar pulse is obtained 




-max 


(2t - t h ) l f 2 


(1-4) 


It differs from (1-3) by the inverse square root time dependence of the 
velocity. 

Another and simpler breakdown criterion was proposed in 1965 by 
Ambartsumyan and al. [44]: if breakdown occurs at the lens focus for a 
laser power P n , it will occur at station r for a power P equal to P r 
times the ratio of the beam cross sectional area at r to that at the focus. 
Again for a triangular pulse the wave velocity comes out to be 
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D = K' (isssy.) 

c max 


1/2 


T/2 


(1-5) 


tn a t 

where K' is a constant containing P n as well as the focal cross section 
area. Canto, Reuss and Veyrie in 1968 [45] introduced the cascade time t 
in a similar model: breakdown occurs at station r and time t if the heat 
flux at r and time t - t was a critical value J g . The wave velocity is 
modified from (1-5) only by changing t into ( t -t) They also 

pointed out the possibility for transparent gases of a breakdown wave 
based on the transmitted power traveling in the direction opposite to the 
lens. 


In summary, as shown by (1-3) to (1-5), whatever the breakdown 

criterion used, the wave velocity is large for steep pulses (large — ) 

tmax 

and long focal lengths (small tna) ; the time dependence is more 
sensitive to a detailed description of the breakdown mechanism. 

The radiation transport wave was also proposed by Raizer in 1965 [37]. 
This elegant idea has not been subject to further developments since then. 
Raizer observed that in air, although the hot plasma is transparent to the 
thermal photons it emits, these have a very short free path in the cold 
gas ahead of the front. Therefore an ionizing precursor is formed; when 
the level of ionization in this precursor becomes sufficiently high the 
laser light is intensely absorbed and a new layer of hot plasma forms. 

The whole process is continuous and can be characterized as a radiation 
diffusion process involving two radiation absorption lengths: the mean 
free path of laser photons in the hot gas and the mean free path of thermal 
photons in the cold gas. Raizer computed the forward emitted thermal 
radiation from the plasma by assuming its shape to be a semi-infinite 
cylinder the cross section of which is the focal cross section of the beam. 
He arrived at a velocity dependence on heat flux 

The three mechanisms for the expansion of the plasma under the 
influence of the laser beam are in competition and it is believed that the 
fastest velocity corresponding to a particular physical situation is the 
one which is observed in an actual experiment. Slow rising pulses and 
short focal lengths are associated with radiation driven detonations whereas 
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fast rising pulses and long focal lengths are explained by a breakdown 
wave mechanism. Indeed in an experiment by Floux, Guyot and Langer in 
1968 [46] a laser pulse was used which exhibited two distinct slopes in 
the rising part; the wave speed corresponding to the steep slope was shown 
to agree with a breakdown wave model whereas the wave speed associated with 
the lesser slope agreed with the detonation theory. The only case where 
the radiation transport wave has been used is by Raizer [37]; he found that 
for the particular physical situation corresponding to the experiment of 
Mandel'shtam [14] the radiation wave velocity came out very close to the 
detonation wave velocity. 

2 . Present Work 

As seen from the above bibliographical survey, theoretical studies of 
the plasma expanding under the action of the laser beam have been dealing 
mostly with the wave front and have been limited to one-dimensional models. 
The only descriptions of the flow inside the plasma have been for a 
spherically symmetric spark whereas in experimental conditions oblong 
asymmetric shapes were reported. From the point of view of nuclear fusion 
the geometrical and temporal distributions of physical variables are 
important, in particular temperature and density. As far as plasma 
containment is concerned an understanding of the spark shape as a function 
of the laser power distribution would be helpful. 

The present work is concerned with the fluid mechanical theory of the 
developing spark during energy deposition from the laser beam; the wave 
front is assumed to be a radiation-driven detonation wave. Consider the 
cylindrically symmetric geometry of fig. 1: the laser beam focused at the 
origin 0 has a certain heat flux distribution in the azimuthal coordinate 
0. A non-spherical plasma is formed and grows with time being led by the 
detonation front. The purpose of this thesis is to investigate the flow 
inside the spark as well as the front shape dependence on the angular 
distribution of the laser power. The general equations describing the flow 
and the boundary conditions corresponding to an energy absorbing narrow 
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layer are presented in Chapter I. A self similar analysis is proposed as 
well as conditions required for it. In Chapter II a perturbation scheme 
is derived for a small perturbation of the laser power from a purely 
spherical distribution. The zeroth order spherically symmetric solution 
is reviewed and the first order perturbation from it is obtained in a 
series of Legendre polynomials in cos 6 . Numerical integration from the 
origin to the boundary was used. The dependence of the detonation front 
shape on the power distribution is presented together with velocity, 
density and pressure profiles inside the spark. This first order perturba- 
tion solution exhibits a singularity at the focus which is analyzed and 
discussed in Chapter III. Finally, Chapter IV contains a summary and some 
concluding remarks. 
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CHAPTER I 


GENERAL EQUATIONS AND BOUNDARY CONDITIONS 
1. Boundary Conditions for an Absorbing Layer 

The plasma boundary is described as a moving energy absorbing layer 
constituted by a high velocity curved shock wave followed by a narrow 
region in which the energy carried by the laser photons is deposited. The 
laser beam is focused on the origin 0 and is characterized by a 6 - distri- 
bution of heat flux; the power may also be time dependent. The strong 
shock elevates the gas temperature and density so that the laser photon 
mean free path becomes very small (Raizer [37] mentions for air at normal 
density and temperature of - 5 10^°K mean free paths of the order of 5 10 ^ 

cm.). Rankine-Hugoniot jump conditions will be derived from the undisturbed 
gas "o" to the hot gas behind the absorption layer "d". Three cases are 
considered: normal steady wave, oblique steady wave and unsteady wave. 

The following assumptions are made: a) the thermal reradiation from 

the plasma is neglected on the basis that the emission length (typically 
several cm for conditions above) is much larger than the plasma dimension 
(a few mm at the end of the pulse) thus making the plasma transparent to 
its own thermal emission; b) the layer width is small compared to the 
front radius of curvature and to any length associated with changes of 
heat flux along the surface: this allows the application of locally one 
dimensional jump conditions c) the perfect gas law is used which is not 
an unreasonable description of a fully ionized plasma d) radiation pressure 
and energy are neglected; Zeldovich and Raizer ([47] p. 141) show that in 
air at normal density, due to multiple ionization, temperatures of 0(3.10^ 
°K) must be reached for the equilibrium radiant energy to become comparable 
to the fluid energy. 


d 


Pd 

u d Pd 
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Steady normal wave : 


Po 

Po 




Consevation of mass, momentum and energy from ”o" to "d" yields the 

following, set of jump conditions: (see e.g. Raizer [37]) 

p u = p' u 
o o d d 

p + p u 2 = P, + P ,u 2 


o o 


d d 


— u + e + P° + 

20 o 

Po ^o“o 


J =1 u + e + Pd 

PnUr, 2 P d 


( 1 ) 

( 2 ) 

(3) 


p is density, u is velocity, p is pressure and e is internal energy per 

unit mass. It is noted that J/p u is the energy per unit mass associated 

oo 

with the laser heat flux J. To obtain boundary conditions on the physical 

variables at "d" equations (1) to (3) are solved for u., p., p, making use 

d d d 

of the perfect gas law and defining the Mach number ahead of the wave 


M 


= p ° u o 
YP_ 


(is the ratio of specific heats at constant pressure and 


volume) 


u = u f 
d o 

r-1 

P . = P f 
d o 

2 

P. * P„ + P u (1 - f) 
a o o o 


(4) 

(5) 

( 6 ) 


from which the Mach number behind the wave can be obtained as 
f 


<■ 


Y[- 


,M 2 

' o 


+ 1 - f] 


where f is defined by 
f = 



1 2 2 -i- 

d-1) " 2(Y " 1J pu 3 

M o o 

o 


(7) 


( 8 ) 


When J = 0 the + sign in formula (8) corresponds to the trivial 
solution whereas the - sign gives the familiar density jump across a shock 
wave. The necessity for the radicand in (8) to be positive yields a 

minimum possible wave velocity: for example for a very strong wave M one 

J 1/3 ° 

gets u > [ 2 (y2 _ 1) — ] . The equal sign corresponds to the Chapman- 

o ~ po 

Jouguet detonation velocity obtained by Ramsden and Savic (1-1). The 

choice of the - sign in (8) gives an overdriven detonation characterized 

by a subsonic downstream flow (M, < 1); the + sign corresponds to an under- 

d 

driven detonation with a supersonic (M, > 1) downstream flow. For constant 

d 
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upstream flow conditions and varying heat flux it is seen that the gas 
velocity behind an overdriven detonation is less than behind a Chapman- 
Jouguet detonation which in turn is less than that behind an underdriven 
detonation; on the contrary density and pressure are larger behind an 
overdriven detonation than behind a Chapman-Jouguet one, the latter being 
larger than those behind an underdriven one. 

In the model adopted here only Chapman-Jouguet or overdriven 
conditions are physically acceptable; behind the leading shock wave the 
flow is subsonic and further energy addition can only increase the Mach 
number towards one. The crossing of the sonic point is impossible unless 
heat is somehow removed in the supersonic region (e.g. by radiation losses). 
It should be noted that for continuous changes from "o" to "d" (i.e. 
without shock) like a breakdown wave or a radiation-supported wave the 
- sign would be the appropriate one and thus would correspond to an under- 
driven deflagration. 

Steady oblique wave : 



The upstream velocity, parallel to the heat flux vector forms an 

angle a with the direction normal to the wave front. The jump conditions 

(1) to (3) apply now the normal velocity and heat flux components and the 

velocity component v parallel to the front is preserved. Boundary 

conditions are obtained, with u , the normal velocity behind the wave: 

nd 

u = u cos a f (9) 

nd o 

v, = u sin a (10) 

d o 

P = P / f (11) 

a o 

p = p + p (u cos a) 2 (1 - f) (12) 

d o o o 
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The total Mach number behind the wave is 
2 


M d = 


c 2 ^ 2 
f + tn a 


( 13 ) 


Y f l 


+ 1 - f] 


M 


on 


M is the upstream Mach number based on the normal velocity component 
on 

u cos o; f is now defined as 
o 


f = 


Y + 1 


[Y + 


M 


on 


- 2 (Y 2 - 1) 

M p u cos a 

on o o 



(14) 


The Chapman-Jouguet wave is again defined by the radical in (14) being 

zero. This can be seen to correspond to a normal downstream Mach number 

2 J 

of one. A minimum possible wave speed is associated u = [2 (v -1) 

o M o 

1 1/3 

°°)which is larger than the corresponding normal 


] 


(for M 


2 

cos a 

minimum wave velocity. The flow behind a Chapman-Jouguet wave is supersonic 

(M > 1) and can be supersonic, sonic or subsonic behind an overdriven wave 
d 

depending on the magnitude of the radical. 

Unsteady wave : 

o 

_ D 


u = o 
o 

P„» 

o o 


a7) 
/ 


J 


/ 


P d’ P d 


The upstream flow from the two steady cases is brought to rest by 

superposing a velocity -u . Regardless of whether the steady wave was 

o 

normal or oblique, the unsteady velocity u, behind the wave is normal to 

d 

the front; D is defined as the normal component of the front velocity. Of 
course the heat flux will generally form an angle a with the normal. The 
following unsteady boundary conditions are obtained, where J represents the 
magnitude of the heat flux: 

(15) 

(16) 


u d = D (1 - f) 


P = P / f 
d o 

P, = P„ + P n V* 
d o o 


(1 - f) 


(17) 
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Y + 1 


[ y + ^7 ±Ja - ~> : 2 - 


M 


M_ 


2 (Y 2 ~ 1) ^f-J 

P o D 


(18) 


M is the Mach number based on D, p and p . It is 
D oo 

noticed that contrary to ordinary shocks f does not go to a constant when 

M -*■ » but generally remains a function of the wave velocity and front 

2 

location (through J -1/r and cos a). This feature makes boundary conditions 
rather more difficult to apply than for ordinary blast waves. 


2. Basic Differential Equations 

The differential equations describing the flow behind the absorption 
wave are those for compressible non-viscous, non-heat conducting fluid. 
Neglecting thermal reemission makes the flow adiabatic and hence isentropic. 
Perfect gas law is assumed to apply. Equations are simultaneously 
presented in spherical and cylindrical coordinates; r is the spherical or 
cylindrical radius, 0 the azimuthal angle measured from the axis of 
symmetry in spherical geometry and a fixed reference direction in cylindri- 
cal geometry; u is the radial velocity component, u the tangential veloc- 
r 0 

ity component. Axial symmetry exists in spherical geometry and physical 
quantities are independent of the z - coordinate in cylindrical geometry. 

No swirl is present. If v is defined such that v = 1 applies to the 
cylindrical case and v = 2 to the spherical case, the continuity, r - 
momentum, 0 - momentum and isentropy equations are: 


3 p 3 v 13 v-1 

^ + 17 (pu r ) + 7 PU r + 7 i? (PU 0 ) + T COt6 Pu 0 


= 0 


9u 9u . 9u 

— J- + u — - + — u a — - 
r 3r r 9 30 


3t 

3 “e x 3u e l 3u e 

iT "tiT r'lF 


_0_ 

r 


1 3£ 

p 3r 


, 9 . 3 . 9 

( + u r i7 + u e i? 


+ U r U 9 _ _ 1 1 i£ 
r p r 30 

)( - £ ) = 0 
PY 


(19) 

( 20 ) 

( 21 ) 

( 22 ) 


3. Self-similar Processes 

The solution of the set of four non-linear partial differential 
equations (19) to (22) in three independent variables r, 0, t subject to 
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boundary conditions (15) to (18) which not only apply at an a priori 

unknown location r (0, t) but also contains the unknown quantity 

d d 

through D, cos a, J, appears a formidable task indeed. Furthermore, 
initial conditions have to be supplied which depend on the complicated 
starting process, namely the initial breakdown mechanism briefly discussed 
in the introduction. Under appropriate conditions self-similarity will be 
shown to apply, this allows to reduce the number of independent variables 
to two: 0 and the similarity variable X Ir/t* 1 . In a self-similar process 
the values of the physical variables describing the flow at station r and 
time t can also be used to describe the flow at earlier time and larger r 
or later time and smaller r, provided the value of 1 remains the same. 
Self-similar flows do not need initial conditions and because of that 
generally fail to properly include the initial stages of a physical 
phenomenon . 

Spherical and cylindrical cases will be treated simultaneously. The 
laser beam is characterized by its angular and temporal power distribution 
(pi 0, t) from which the heat flux hitting the wave front is 

J ( 6)t ) = (23) 

[r d (6,t)] V , v , • . 

Self-similarity is obtained when the two following conditions are 
satisfied: a) p is specialized to a power- law time dependence 

Pi 0 , t ) = P L tf g (0) (24) 

where P is a constant representing power per unit time to the k and g (0) 

is a non-dimensional function of the azimuth; b) very strong waves are 

considered i.e. -*• <*> and p^ r*- 0. This in fact means neglecting the 

static pressure ahead of the wave compared to its dynamic pressure; for high 

velocity waves as those experimentally observed this assumption is 

certainly valid. An analysis is made in Appendix A which includes finite 

pressure p for a spherically symmetric spark: the results are presented 
° 2 

in the form of a perturbation in 1 / tl from the self-similar solution. 

D 

Dimensional arguments can be used to determine the proper similarity 
variable (see Sedov [42] ). The overall physical problem possesses only 
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two dimensional parameters: P of dimension mass times length to the v 

L 

divided by time to the (k + 3) and the undisturbed density p of dimension 

o 

mass over length cubed. It is seen that the ratio P, / p contains only the 

L o 

dimensions of length and time. The non-dimensional similarity variable is 
thus defined 


< + 3 . . r P. 

n = — p— and A = [ _L (y-l)aj 

P 

o 


1 

v + 3 


a is a constant to be obtained as part of the solution. It will be chosen 

such that the boundary of the spherically symmetric spark is X = 1. From 

(25) it is seen that the boundary is now located at 

X = h (0) (27) 

or in physical variables r, = h (0) A t n . 

d 

The self-similar variable can also be obtained by considering the 
boundary conditions; in particular it is necessary that the density ratio 


across the wave be time- independent . Looking at the lastt^rm of the radi- 
cand of (18) it is required that the quantity — (the dot 


v • 3 

r. (r,) J 
d d 


represents time-derivative) be not a function of time. Integrating that 
^ ... . . ic + 3 


quantity with respect to time one obtains r . 


t — — - ; defining X as r/r, 
v + 3 d 


one recovers (25) - (26) when all constants are taken into account. 

The explicit time dependence in the set of differential equations (19) 

- (22) can be shown to cancel out when the velocities, density and pressure 

are respectively non-dimensionalized by A t" \ p and p (A t* 1 ) ^ (see 

o o 

e.g. Rae and Kirchner [48])* This result actually holds for any values of 
A and n. This type of self-similarity although limited to spherically 
symmetric geometry was first proposed by Taylor [49] for ordinary blast 
waves, it differs slightly from Sedov's [42] type in which the non- 

2 

dimensionalization is carried with respect to r/t, p and p (r/t) . 

o o 

The self-similar equations are obtained by changing variables in 
equations (19) — (22) from r, 0, t to X,0 and introducing the self- 
similar physical quantities (represented by capital letters) functions only 
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of A,0 : 

= n A t n V ( A,0) : 

u =nAt n ^W ( A , 8 ) 

0 

p = p R ( A,0 ) * 

O 4 . j 

,p = p (n A t n ^") 2 P ( A,0 ) 
o 

canceling out the explicit t - dependence in all equations one gets 

D « R + R [V. + v 7 + ~ W + (v -1) cot0 7 W] = 0 
Ac A ^ A A 8 A 

D 'V+^-^V- — +- P =0 
A0 n A R A . , 


A0 


W + — w + ^ + ± i P 
n A R A 0 


(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


°A0 * p + 2 P + Y P [ V A +v f + I w e + ( V ~D cote J W] = 0 (35) 


where D 


AS 


D. 


A 0 


is a linear differential operator defined by 


(36) 


The boundary conditions (15) - (18) are also recast in self-similar form. 
First, from (27) the normal component of the wave velocity and the angle a, 
defined positive from the external normal direction to the front to radial 


direction, turn out to be: 

1 , »2 -'1/i 

D=nAt nl (l+^ L ~) 

, 2 h 

cos a = ( 1 + - — J -1 ^ 2 

h 2 2 

sin a - -l (1 + ^) 1/2 

n 


(37) 

(38) 

(39) 


h' is the derivative of h ( 8 ) with respect to 0 . Then, remarking that the 
velocity behind the wave which is perpendicular to the front yields both a 
radial and a tangential component the final boundary conditions come out: 
at A . h (9) 


) 


-1 


V - (1 - f) h (1 + - _ 
d h 2 

w = -(1 - f)h’ (1 + ~) 1 
d h 2 

R d ■ 7 ' , 2 

P = (1 - f) h 2 (1 + ~ - ) 1 
d h 2 


(40) 

(41) 

(42) 

(43) 
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Only the - sign will be retained in f, consistent with the chosen detonation 
wave model. The problem has now been reduced to a system of four non- 
linear partial differential equations in two independent variables. The 
boundary conditions still contain an unknown function h (8) . 

It may be noticed that the self-similarity presented here is the 
same as for spherically symmetric geometry; this results from the fact that 
the angular distribution of power g (8) does not introduce any new physical 
length in the problem. 
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CHAPTER II 


SMALL PERTURBATION OF A SPHERICAL SPARK 
1. Small Perturbation Assumption 

Throughout this chapter only constant power in time will be 
considered i.e. k = 0 and hence n = 3/5. The geometry will be spherical: 
v = 2. Consider the following arrangement: a battery of lasers all 
focused at the same point creates a slightly non-spherical angular 
distribution of power: 

g(0) - ^ I 1 + eG (6) ] (45) 

where e is a non-dimensional parameter much less than one. The object 
of this investigation is to determine the spark shape and the plasma flow 
resulting from such a power distribution. Since a uniform change in power 
level can easily be treated by a perturbation of the definition of A 
(see eq. (25)) the function G (0) will be chosen such that the total power 
P L remains unchanged. This is expressed by the condition 

F G (0) sin 0 d0 = o (46) 

It is assumed that the displacement of the boundary of the 
initially spherical spark, resulting from the imposed power perturbation 
is also proportional to e ; the wave front is thus 

A - h (0) = 1 +e H (0) (47) 

where H (0) is an a priori unknown function to be found as part of the 
solution. 

In the light of equations (45) and (47) the boundary conditions 
previously derived (40) to (44) can be expanded in a power series in e . 
First it is observed that, provided the radicand in (44) is different from 
zero in the spherically symmetric case, and this will be shown later in 
this chapter, the function f can be expanded as 

f = f (o) + e K [G (0) - 5H (9 ) ] + 0(e 2 ) (48) 

where f^°^ and K are constants depending only on the spherically symmetric 
solution (which from now on will also be referred to as "zeroth - order 
solution") defined by 
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(49) 


(o) 


— [y - (1 - ^-> ' ] 

Y+l 2irn 2 3 a 


y +1 " 1/2 

and K « 7Z~i (1 - 3 — ) (50) 

47Tn a 2Trn a 

Expanding in powers of e equations (40) to (44) yields the boundary 

O 

conditions where terms of 0 (e^) have been dropped: at X = 1 +e H(6) 
V d = 1 - f<°) + e [-K G(9) + (1 - f<°> + 5K) H(6) ] (51) 


W, 


-£ (i - f(°>) H'(e) 

1 K 


f(o) 


- e 


.(o) 


(G(0) -5 H(0) ] 


(52) 

(53) 


P d = 1 - f (o) + e[ -K G(6) + { 2 (1 - f (o) ) + 5 K> H(0) ] (54) 

The form of the boundary conditions (51) to (54) naturally suggests 
expansions in powers of e of the physical variables everywhere in the flow 
field 

r~\ fi \ 

(55) 


V (A, 6) = v (o) (A) + e V (1) (A,0) 


W (A,0) = e W 


( 1 ) 


(A,S ) 


R (A,0) = R<°> (A) + e R (1 > (A,6) 

P (A ,6) = p(°) (A) + e p(l> (A, 8) 


(56) 

(57) 

(58) 


The zeroth-order functions representing the spherically symmetric case are 

dependent only on A . It is noticed that, although the expansions of the 

boundary conditions are exact, provided that assumption (47) is valid, 

there is no guarantee that the expansions of the field variables (55) to 

(58) are uniformly valid. In fact it will be shown that they do break down 

in the neighborhood of the focus: A = 0. 

Proceeding with the substitution of (55) to (58) into the flow 

equations (32) to (36) one collects the 0 (1) and 0(e) terms to obtain 

the following two systems of differential equations: 

(V(°) -A) r(°)’ + R<°) V<°>’ + 2 R - 0 (59) 

n^l v (°) + ( V (°) _ A ) v (o)’ + 1 p(°)’ = o (60) 

n r(°) 

2 Bzi p(°) + (y(°> - A) P (o) ’ + yP^ + 2y P ~— = 0 (61) 

n A 

In the above zeroth-order non-linear system of ordinary differential 
equations * stands for — — - . 
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The first order system is a linear system of partial differential 
equations, the non-constant coefficients of which are function of the, 
zeroth-order variables. It reads: 


R<°> V™ + (R (0 >’ + 2 ^ ) V (1) + ± R (0) (W^ + cote W (1) ) + 
(V (0) - A) R^ + (V (0) ' + 2 j ) R (1) = 0 

(v (o) - a) v ( 'J ) + (— + v (o) ) v (1) - R (1) + ~~ p ( 0 

A n _(o)2 „(o) A 


< v <°> . „ „W) + <8=1 + 

A n A 


(o) 


) w (1) + 1 


, , 1 P (1) -0 

(o) a . , e 


( 62 ) 

(63) 


(64) 


y P (o) V < } > + (P 


(o)’ .. - P ( °\ „(1) 1 (o) ,(1) ^ . (1). . , 

x T *.*• + 2 y — ) v +Y y p (W e + cot- e w ) + 

(V (o) -A) P ( ^ + (2— +Y V (o) + 2Y ) P (1) = 0 (65) 

An A 


The boundary conditions on the two systems are readily obtained from 
(51) to (54); for convenience, they will be transferred to the fixed loca- 
tion A = 1 by means of Taylor expansions around A = 1 +e H (0). One thus 
obtains the following zeroth order and first order sets of boundary 
conditions: 

(o') ( o'] 

, ...... ( 66 ) 

(67) 


v (o) (i) = l-f (0) 


R (o) (1) - 1 


P (o) (1) 


(o) 


1 - f 


(o) 


( 68 ) 


and 


v (1) 

( 1 , 

0 

) 

= - KG( 0) 

+ 

[1 

- f ( °> 

+ 5K 

- v (o) ’ (1)] 

H(0) 

(69) 

w (1) 

( 1 , 

0 

) 

- - U-f Co> ) 

H’ 

(9) 




(70) 

R (1) 

( 1 , 

0 

) 

K 

9 

G 

(0) 

+ [5 ■ 

K 

9 

■ R (oV (1)] 

H (0) 

(71) 





f <°V 




f to) . 




P u> 

(1, 

, 0 

) 

= - KG (6) 

+ 

[2 

(1 - f 

( °>> + 

5 K -,P (o) 

(1) 3 H- (0) 

(72) 


2. Spherical Spark - Zeroth-order Solution 

Although the spherical spark problem was treated by Champetier and al. 
[40] and Wilson and Turcotte [41a] a new investigation was prompted by the 
necessity of a readily available solution in order to compute the coeffi- 
cients of the first order set. It was also observed that a slight 
discrepancy exists between the pressure profiles obtained by the above 
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investigators: Wilson's and Turcotte exhibit a dip around X - 0.5 which is 
not present in Champetier's and furthermore Wilson and Turcotte's pressure 
at the origin is about fifteen percent higher than Champetier's. The 
density profile in Champetier and al. [40] appears to be erroneous; al- 
though they correctly indicate that the wave front is overdriven, they 

mention a density ratio across the wave: f = y which, of course, is 

y+ 1 

that of a Chapman- Jouguet detonation. 

The approach used here is quite different from the previous investiga- 
tions: the set of equations (59) - (61) is directly integrated as a 
function of A from the origin out to the boundary. This necessitates a 
detailed study of the behavior of the physical variables in the neighborhood 
of the origin; the understanding of this singularity at the origin also 
throws some light on the more complicated singularity arising in the first 
order perturbation problem. . Since the boundary conditions (66) to (68) 
with f^°^ defined in (49) contain the unknown constant a, and hence 
provide only two known boundary conditions, expansions near the origin are 
needed which contain only two independent constants. 

I 

2,1 Expansions near A = 0 

For convenience a change of variables is made: 


V> = x ° (73) 

vf 0) = A y 1 (u) . (74) 

R (0> = ~ y, (p) (75) 

A 2 z 

P (o) = y 3 (p) (76) 

where a is, at this point, an arbitrary constant. Transforming equations 
(59) to (61) and canceling out all explicit A - dependence, one obtains: 
(with ' = -r ) 

dp ; 

r, (y x + 2) y 2 + ap[y 2 y^ + (y^-l) y^l = 0 (77) 

■ y x (y x - y' 2 + a d[(y 2 - D y 2 y^ + yj] = o (78) 

(3 Y y 1 + 2 ^) y 3 + a p[ y y^ + ^ - 1) y^] = 0 (79) 

Expansions near A = 0 are sought of the form: 

1 2 3 ’ 

y 1 = a + a p + a p + a p + (80) 
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(81) 


2 3 

y 2 = b l P + b 2 P + b 3 P + “ * 

2 3 

y, = c + c u + c w + c u + 
3 o 1 2 3 


(82) 


Note that the constant term in the expansion of y being zero renders the 
density It bounded at the origin. The leading term in equations (77) to 


(79) yields respectively: 


[a + 2 + a (a - 1)] b = 0 
o^o 1 

a (a - — ) b + a c = O' 
o o n 1 1 

(3 Y a +2 — ) c = 0 
o no 


(83) 

(84) 

(85) 


Therefore the expansions near the origin contain two arbitrary constants 

( 86 ) 


b and c„ , the value of a is found and a is determined 
o 1 o 

= _ __2_ n-1 _ _ a + 2 

a o 3 Y n ~ ” ° 


a -1 
o 


c is readily obtained from (84): 

/ 1 \ 1 , 
c. = - (a - — ) a — b, 

1 o n o a 1 

Numerical values corresponding to Y = 5/3 (and recalling n = 3/5) are 


(87) 


a^ = 4/15 and a = 34/11. a being non-integer makes all functions singular 

at the origin in the sense that although the function itself is bounded, 

(o) M 

there always exists some derivative which is not (e.g. K is infinite.) 

The self-similar velocity, density and pressure are seen to have a 
leading term respectively proportional to A , ^2/11 an{ j a constant. From 

the definitions (28) to (31) the r, t dependence of the physical velocity, 

. . . _ . . J , . 12/11 - 12/11 
density and pressure are obtained respectively as r/t, r . t , 

-4/5 

t . This kind of flow where the velocity is proportional to r and the 

pressure is independent of r is usually referred to as uniform spherical 

expansion. It is worth noting that the temperature is unbounded at the 

origin ( r . t^^"*) as it indeed must in a self-similar spherical ' 

non-heat conducting process. The initial stage of a spark is characterized 

-2/5 

by an infinite front velocity, D ~ t ; this raises the entropy of the 
gas particle at the focus to infinity. But, by symmetry this particle stays 
at the focus so that entropy cannot be convected away from the origin; on 
the other hand it cannot diffuse away since non-heat conducting fluid is 
assumed. Thus the temperature at the origin has to remain infinite at all 
times . 


22 



In order to improve the accuracy of the numberical solution several 
more terms are computed in expansions (80) to (82). For Y = 5/3 the 
following numerical coefficients are obtained where B = b^/c Q : 

y = ~ + 2.69692 10~ 2 8 V + 3.10402 l(f 3 B 2 V + 2.01178 10~ 4 B 3 p 3 (88) 

y = b u (1 + 8.54504 10" 2 6 yt- 1.2256 10~ 2 B 2 V + 1.56620 lo" 3 
2 1 3 3 

B V ) (89) 


y„= c (1 + 1.20784 10 _1 6 v + 1.99936 10 2 
3 o 

„3 3 . 

B v ) 


8 2 v 2 + 2.91446 10~ 3 


(90) 


These expansions containing two arbitrary constants are sufficient to 
start a numerical solution which has to satisfy two boundary conditions. 
However, it could be argued that the system (59) - (61) is third order and 
a third solution has been implicitly discarded. In fact, by inspection a 
solution depending on one arbitrary constant C can be found (for Y = 5/3) 


v (°) = 1 x 
6 


R (0) = C X 15 


» 


25 _1 
36 17 


C> X 


17 


(91) 


This solution has to be discarded indeed since it does not have the required 
physical feature of an infinite temperature at the origin. 


2.2. Numerical Integration and Results 

The system of three non-linear ordinary differential equations (77) to 

(79) is integrated from 0.01 to 1 using a standard fourth-order Runge-Kutta 

' ’ /' > - *’ • “ _2 

technique with a step-size of 10 or less. Expansions (88) to (90) are 

v . ■ - * 

used to move away from the origin. Two boundary conditions are to be 
satisfied, for example: 


y i ^ = y 3 ^ ; y i ^ = 


y 2 ..U) 


(92) 


The values of the constants b, and c are obtained by a linear correc- 

1 o 

tion procedure: considering the y's to be function not only of y but also of 

N N 

b and c , their values corresponding to a new guess b , c are related to 
1 ; o - v ±o 

o o 

the values corresponding to the old guess b. , c by a Taylor expansion: 

3 1 ° 


y (y; b^, c^) = y (y; b°, c°) + 3 - ? - (y ; b° c°) H 

1 o 1 o ° b, 1 o 13c 


/ l° o, . 

( M ; b , c ) A c 

1 o o 

No No 

with 4b, = b, - b and Ac = c - c . 

Ill o o o 


(93) 


The system (77) to (79) and the 
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expansions (88) to (90) must be supplied with six extra equations and six 

extra expansions for the derivatives of the y's with respect to b, and c ; 

1 o 

these are obtained simply by taking derivatives of (77) - (79) with 
respect to b^ and c q as well as of equations (88) to (90). The full system 

j O ■ 

of nine equations is then integrated for an initial guess of the pair b^, 

c°. The values of the y's and their derivatives for the initial guess are 

now known at the boundary. Application of the boundary conditions (92) on 

the values of y's at the new guess gives a system of two inhomogeneous 

algebraic equations for A b^ and A c q . Solving this the A's are found and 

hence the new guess. The procedure is repeated until both a's become small- 

-12 

er than a prescribed value, namely 10 in the present case. This 

correction scheme was found to converge fairly rapidly requiring fifteen 

iterations to bring b and c from an initial guess of 1.0 to their final 

1 o 

values : 

b = 0.9299591 c = 0.173074414 

1 o 

The self-similar velocity, density and pressure profiles obtained are 
plotted on figure 2 as a function of the self-similar variable X. The 
boundary values are 

R (0> (1) = 1.85043402; V (o) (1) = P (o) (1) = 0.459586242 
where all digits are significant. From that the constant a related to the 
wave front speed can be found 

a = 1.990 

2.3. Discussion and Conclusions 

A highly accurate numerical solution has been obtained; this results 

from the use of several terms expansions near the origin as well as a severe 

-4 

convergence criterion combined with step-size as small as 5 10 in the last 
iterations. No pressure dip similar to Wilson and Turcotte's is observed 
and the ratio of the pressure at the origin to that just behind the absorp- 
tion wave is 0.37659 compared to Wilson and Turcotte's value of 0.441. 
Aside of these differences the velocity, density and pressure profiles and 
in particular the values at the boundary, agree fairly closely with Wilson 
and Turcotte's results. 

(o) Y 

The density ratio at the front f = 0.540413758 is less than ' 
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thereby indicating an overdriven detonation front. This result is 
important for the perturbation scheme since it justifies the previously 
cited eq. (48). Also note that the whole flow inside the spark is 
subsonic. 


Expansions 


A comparison with ordinary constant energy blast waves is enlighten- 
ing. From a mathematical point of view, their treatment would be very 
similar. Changes of h to 2/5 and of the definition of A are needed and 

the boundary conditions are modified by having f^ 0 ^ — — . Expansions 

y +1 

near the origin are of the same form with a Q = 3/5 and a = 13/2. It is 
observed that the boundary conditions do not contain the unknown wave 
velocity "a" and therefore appear to constitute three independent boundary 
conditions in contradiction with the expansion near A = 0 containing only 
two independent constants. In fact by manipulating the equations (59) - 


(61) the following integral relation can be obtained: 

<V<°> [ P<°> ♦ (Y-l) S <0) » ( ° )2 || X-1 + (y-D p<o> v <°> 


+ *4 ill p ( o) ♦ (Y-l) 2 


1 |x-l + p( ° >v |x.l 

/ \ n 


dA = 0 


It states that the total energy swept in by the moving front plus the 
pressure work balance the change in total energy of the flow inside the 
spark. For a constant energy blast wave (n = 2/5) the integral in (94) 


drops out and hence a relation is obtained solely from the differential 


equations, between the values of the physical quantities at the boundary. 
Clearly only two boundary conditions need be applied. 


Physically the constant power laser spark has a faster moving front 
than a constant energy blast wave which is natural in view of the contin- 
uously added energy. The mechanism pushing the leading shock wave is quite 
different: in a blast wave the energy released initially at the origin is 
redistributed as the blast grows whereas in the laser spark the energy is 
deposited right behind the shock driving it forward. This is the reason 


why the boundary conditions appear different: dependent on velocity for 
a laser spark and independent for a blast wave. The flow near the origin 
is very much alike in the two cases: constant pressure and v-r/t. 
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It should be noticed that for the same y= 5/3 the density profile of 
a blast wave drops much more steeply near the boundary. This, actually, is 
only a manifestation of a very important difference between the two flows, 
namely the fact that when Y~*l a blast wave exhibits a Newtonian layer of 
concentrated mass near the shock, contrary to the laser spark which does 
not. Figure 3 shows the velocity, density and pressure profile of a 
constant power spherical spark for y= 1; they clearly do not have a , layer 
of concentrated mass near the boundary. The procedure to obtain these is 
the same as for y- 5/3 although convergence appeared somewhat slower. The 
existence of a Newtonian layer for blast waves permits the considerable 
simplification of treating the flow as locally radial (i^e. dropping of 9 - 
derivatives compared to r - derivatives) , this possibility was exploited 
by Laumbach and Probsteih [50,51] in their study of blast waves in non- 
homogeneous atmosphere. However, the author disagrees with their use of 
the same simplifying assumption in their treatment of blast waves in in- 
homogeneous atmosphere including the thermal radiation emanating from the 
inside of the blast: the boundary conditions are then modified. to look like 
those derived in Chapter I (equations (A) to (8)) and this should destroy 
the possibility of a Newtonian layer. „• 

The reason why a Newtonian layer is possible or not is readily seen v 

from equations (15) to (18.) : when Ytl». in the case of a shock- wave the 

density tends to infinity at the same time as the velocity u^ tends to D 

thereby trapping most of the mass inside the-, blast wave in a layer next to 

the shock. Therefore for a blast wave with Y^l the velocity profile tends 

(o) , 

to a straight line V = n A which implies that particle trajectories are 
also similarity lines: this feature allows the existence of a constant 
pressure, zero density core inside the blast wave into which the mass flow 
is zero; For a constant power laser spark none of these characteristic 
features hold and hence a Newtonian layer is ruled out. 

In conclusion, the non-existence of a Newtonian layer for constant 
power excludes the possibility of using the locally radial approximation 
and a small perturbation scheme has been proposed instead. 
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3. Non-spherical Spark - First Order solution 
3.1. Fourier Decomposition 

Contrary to the zeroth order set, the first order set of differential 
equations (62) to (65) is linear; this feature is exploited by expanding 
the dependent variables in Fourier series in the coordinate 9. First the 
imposed power addition perturbation G (0) is decomposed in a. series of 
Legendre polynomials in cos 0 : ' 

CO 

G(0) = I Pfc (cose ) A . (95) 

l = 1 

where is the Legendre polynomial of order 1 and the are known 
constants. Note that the requirement (46) that the total power should 
remain unchanged is satisfied since the constant term P q is omitted in 
the series (95) and for & >1, all Legendre polynomials satisfy 

/ P (cos©) sin 0 d 0 * 0 (96) 

0 

The first term of the series P^ = cos9 corresponds to a maximum asymmetric 
perturbation: if A^ >0 the laser beam is altered such that the power is 
Increased by a certain amount at 0 =■ 0 and decreased by the same amount at 
0 = it . Higher order polynomials create a more evenly distributed power 

perturbation. When the order & is even the distribution is symmetric with 
respect to the axis 0 = — in addition to the original symmetry about 0 0,n 

and the whole flow presumably possesses two axis of symmetry. 

The perturbation of the spark shape H (0) is similarly decomposed in: 

H(6) = Z P (cos0) X. (97) 

where the are unknown constants since the shape is to be determined as 
part of the solution. 

In order to achieve Fourier decomposition of the differential 
equations and boundary conditions the following expansions are made for the 
velocities, density and pressure perturbations: 

V (1) (A,0) = Z p (cos9) V* (A) (98) 

1 =i 

W (1) (A,e) = £ — (cose)] W £ (A) (99) 

£=1 
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(100) ' 


( X, 0) = z P n 
«,= 1 1 

(cos 0) R^ 

(X) 

( X, 0) = Z p. 

(cose) p 

(X) 

II 

l 



t (1) 

( 1 ) 

- • - - • (101) 

where the quantities — are solely function of A. Again the constant terms 
P q = 1 have been omitted consistently with jthe remark made previously that 
they represent a purely spherical flow which can easily be treated by a 

t 

perturbation of the definition of X. The expansions (98) to (101) are 
inserted into equations (62) to (65). Making use of the defining differen- 
tial relation for Legendre polynomials , 


d 2 P* 

d02 


+ cot 0 


dP i 
de 


+ l(i. + 1) P 


( 102 ) 


to express the quantity (W^ + cot© W^) as function of P fc rather than 

its derivatives it is seen that P (cos0) can be factored out of equations 

x> , 

(62) , (63) and (65) and that Pn (cos0 ) can be factored out of (64) 

d0 . .! 

which is then integrated with respect to 0 producing an arbitrary constant 
C. Each equation, which is now of the form P (cos 0). function of A, is 
multiplied by P^ (cos0) sin0 (with k = 1,2,3** ■) and integrated from 0 to it. 
Using the orthogonality properties of Legendre polynomials: 


f Q P^ (cos0) . P^ (cos0) sin 0d0 = 0 k^Jl 


(103) 


2k + 1 

and remarking that the contribution of the arbitrary constant C is zero, for 
any value k > 1 one obtains : - 

R (o) vj + (R (0) + 2 ^^v k - ~ (0) k (k +1) W k + (V (o) ~A> \ + 


(o)' + 2 V^) . o 


(V 

A »• , 

(V (0) -A) \ + (*^ + V (0) ) \ - P 


(104) 


^ \ + ifey = 0 (105 > 


(V 


(o) 


t 


i-l 


V<°> 


•X) w k + (rrr + ^— ) w k + 


4-T - P = 0 . (106) 

,(°) X k 


(o) _» ( 0 ) 

YP V + (P 
k 


,(o) 


(o) 


+ 2 Y' 


>v-V 


Yk (k + 1) W k + 


(V^ -x) P k + (2 




v (o) - 

+ * V > P k = 0 


(107) 
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In the above equations ' refers to ~ • The boundary conditions are 
treated likewise: the expansions (95) for G (0), (97) for H (0) and (98) 

t 1 ' 

to (101) for the physical variables are substituted into equations (69) 

to (72); after multiplication by P- (cos9) sin 6 and integration over 0 

k 

from -0 to ".ir one obtains the set' of boundary conditions: 

V (1). - KAj+ f 1 -‘f (o) + 5K'~ V (o) ’ (1)] 3^' (108) 1 

W k (l) = - (1- f (o) ) \ , (109) 

- - Tbr \ + [5 £ -fe)2 - K (o) "(i> ) Xfc (HO) 

P k (l) = - K A + [2 (1 - f (o V+ 5 K - P (o) ' CD ] \ (111) 

The differential equations (104) to (197) and boundary conditions 
(108) to (111) apply separately for each value of k = 1,2, 3... The problem 
has been completely split into its different harmonic components: the 
influence of each constituent harmonic of the power perturbation G (0) 

i • . ; i • ' ■ • i ' • ' ’ ' 

can be investigated independently. The simplification is the result of 
adopting a linear perturbation scheme in e and it is remarked that the 
full set of non-linear equations (32) to (35) is not amenable to such a 
treatment. 

The problem now is. reduced to that of obtaining the solution of the 
system of four linear homogeneous ordinary differential equations with 
variable (and known only numerically) coefficients subject to four non- 
homogeneous boundary conditions . These contain the unknown constant Xj^ 
which, in fact, reduces the number of known boundary conditions down to 
three, by eliminating 3^. The system being fourth order, this implies 
that a solution will have to be discarded on physical grounds. 


3.2. Expansions near A = 0 

The set of equations (104) to (107) will be integrated outward from 
the origin to the boundary. Starting such a solution necessitates an 
investigation of the flow in the neighborhood of A = 0. Only the leading 
term of. the expansions of V , W , R and P is sought; to that effect 

K. K K i iC 
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equations (104) to (107) are specialized to the neighborhood of the origin 

by retaining only the leading term in the expansion of their coefficients. 

Use is made of the previously obtained . expansions of the zeroth order 

physical variables (80) to (82): 

b * a_2 V,' + b a X a ^v - k(k+l) b \ a ~ 3 W, + (a -1) \ 5' + 3 a 
lkl k 1 koko 


\ = 0 c „ (112) 

(a -1) X V’ + (— + a ) V - ~ - X ““ +3 i + ± X +2 p' = 0 (113) 

o k n ok, 2 k. b, k 

b i 1 

(a -1) X W’ + (— + a ) W , + X _a +1 P = 0 (114) 

o k n o k b. k 

c 1 c _ _ 

yc V,' + 2 y t 2 - V, - Y k (k + 1) — 2 W + (a -1) X P ' + Y a 

ok A k Ako k 1 

(a + 3)X a P, =0 (115) 

k 


It is observed that the coefficients of P' and P in equation (115) are 

2 K. K 

respectively multiplied by X and X compared to those in equations (113) 
and (114). These higher order terms may be dropped and equation (115) 
reduces to its first three terms. Equations (112) to (115) constitute a 
fourth order system, hence four solutions are sought, each of the form: 

(116) 

X “k (117) 

X \ + 0-3 (118) 

(119) 


\ = A lk X ^ 


\ = A 2k 


\ 


= A 


3k 


P . A„ A \ + “- 1 
k 4k 


When these are substituted into the set (112) to (115) all terms in each 
equation come out to be of the same order in X . Cancelling out all X 
dependence a system of four algebraic equations is obtained for the con- 
stants A., , A„, , A„, , A., . Since the original differential equations are 
Ik 2k 3k 4k 

homogeneous the algebraic system is homogeneous too and m^ must be root 
of an indicial equation, namely the determinant of the system equals zero. 
After some manipulations in which use is made of the definition of a 
(86) it reads: 
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a-2 


0 


0 


(V 1 ) 

(a o -i) \ “ + a o 0 a 0 (a Q ~ \ + ^ 

0 (a -1) m. + ^ - + a 0 1 

° k n o 

m+2 - k (k + 1) 0 0 


= 0 

( 120 ) 


Developing this determinant yields a fourth order algebraic equation in 
n^. Its coefficients do not contain the constants and c Q entering the 
expansions of the zeroth order quantities so that the values of n^s are 
dependent only on the characteristic quantities a Q and a of the zeroth 
order numerical solution. However, the m^s are seen to be function of y 
and of the harmonic considered k. The fourth order equation has been 
solved numerically for y => 5/3, the results are presented in table 1 for 
the first ten harmonics. A definite pattern for the numerical values of 
the roots seems to emerge for increasing k i.e. for perturbations tending 
to be more and more spherically symmetric: the first root is positive and 
increasing, the last root is negative and decreasing and the pair of 
complex conjugate roots has a real part negative and increasing for the 
first three harmonics and a real part positive and increasing from harmonic 
four on. An important feature of the complex roots is that their real 
part always remains less than one, even for large values of k. This will 
be demonstrated in section 3.4. 


The system of four algebraic homogeneous equations for A^ to A^ k 
is then solved for each root. One obtains, with the superscript i referring 
to each root i = 1, 2, 3, 4: 


k (k+1) lk 

i a -2 

3k , i . . , 

(m k " 1 >< a o- 1 ) 


b i 4 


me + 2 i n-1 i 

toTTT) « a o-i> \ a 11c 


(121) 


( 122 ) 


(123) 


For each root, A^ k has been arbitrarily chosen as the free constant. Note 
that for the pair of complex conjugate roots all A's are complex numbers. 
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harmonic k 

roots: m^ 

1 

1.114 

- 0.6253 +0.5500 1 

- 3.500 

2 

1.399 

- 0.2491 ±0.7055 i 

- 4.537 

3 

1.941 

- 0.0137 ±0.6691 i 

- 5.550 

4 

2.719 

+ 0.0998 ±0.5965 i 

- 6.555 

5 

3.619 

+ 0.1510 ±0.5413 i 

' 

- 7.557 

6 

4.569 

+ 0.17.67 ±0.5047 i 

- 8.558 

7 

5.540 

+ 0.1912 ±0.4802 i 

- 9.559 

8 

6.522 

+ 0.2002 ±0.4632 i 

- 10.56 

9 

7.510 

+ 0.2061 ±0.4511 i 

- 11.56 

10 

8.510 

+ 0.2103. - 0. 4421 i 

- 12.56 


Table 1: Roots of the indicial equation for harmonics 1 to 10 

and Y = 5/3. 

From this, the expansions of V^, W^, R^, near A = 0 can be 
obtained as a sum of four linearly independent solutions. For 
convenience the notation is simplified to: m^ refers to the positive root 
and is associated with the real constants A£s; m£ refers to the negative 
root and is associated with the real constants D^s; m rk ± i m ^ k are the 
complex conjugate roots associated with the complex constants B's. The 

, K 

expansions of the k c harmonic read: , 

\ <» 

\ < X) 

\ (X > 

? k <»> 


= A. 


Ik' 


A "k + B lk * + B* k x ^ ik + Plfc A > • <«4) 

m — m 

m-v+im., * 

= A,. A ~ + B„. A lk + B„, A 


®k 

'2k A ‘ "2k 


rk ik 


m 


2k 


+ D 2k A k 


(125) 


= A,. A m k + °- 3 + B„. A ® r k +a " 3+im ik + B* A ra rk +n - 3 " im ik 


'3k 3k 

+ °3k X "k + 

* *4k ' " k+ °^ 1 + B 4k» “ tk 
+ D 4k X K + »-‘ 


. 3k 


+a -l+im. k +B * m rk +^l-l m . k 
4k 


(126) 


(127) 
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* stands for complex conjugate. Note that in each relation the second and 

third terms always combine to form a real quantity as they should since 

they represent real physical quantities. The free constants are now A,, 

lk’ 

Bfk , and the A^ k , B^, Dj k (j = 2,3,4) are respectively related to the 

free constants by (121) to (123). The pair of middle terms, containing 
B's can also be rewritten as: 


B X m rk +im ik + B * A™ 1 * imik = 2A mrk ( B r cos B 1 sin £) 


• 5 = m ik log X 


(128) 


(129) 


where B^ and B i are respectively the real and imaginary parts of B. In 
this form the oscillatory character of the expansions is clearly displayed. 


3.3 Numerical Integration and Results 

The system of differential equations (104) to (107) is integrated from 
A = 0.01 to A = 1 using a standard fourth-order Runge-Kutta numerical 
scheme. The three boundary conditions are obtained from (108) to (111) 
by elimination of X^. The solution is started using expansions (124) to 
(127) in which the terms corresponding to the worst singularity (i.e. the 
negative real root m^) are discarded by setting D lk = o and hence D 2R = 

D^ k = D^ k = 0. This will be justified in the next chapter on the physical 
ground that this singularity represents a source of mass and energy at the 
origin. 

The linear character of the differential system and boundary conditions 

is exploited to determine the values of the free constants A^ k , B^ kf , 

®lki ^ application of the boundary conditions. First a solution is 

obtained starting with expansions (124) to (127) specialized to A^ k = 1, 

B^ kr = 0, = 0. Call this solution where is a column 

matrix V k> w k » R^, P . A second solution S^^) i s obtained with A^ k = 0, 

B =1, B„, = 0 as well as a third solution S, (3) f or a , = 0, B = 0, 
lkr . lki k lk lkr 


B 


lki 


= 1. The full solution is then sought in the form of a linear 


combination of S , S ^) # s : 

k 1: k 


Sk = A lk S R a) + B, , — S 


lkr "k (2) + B lki S k (3) 


(130) 
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The value of the constants are obtained by applying the boundary conditions 
on the full solution S^. Since these are non-homogeneous, a third order 
non-homogeneous linear algebraic system is to be satisfied by A , B and 

IK JLKX 

B,, , and the velocities, density and pressure profiles are obtained as 
lki 

(1) , S k < 2 >, S k (3). As 
a check of the numerical accuracy the equations were integrated over again 
with the correct values of the free constants and it was found that the 
boundary conditions were accurately satisfied. 

Note that the system (104) to (107) was supplied with the set of 
zeroth order equations (59) to (61) so as to have the numerical values of 
the coefficients at all grid points. The seven equations were integrated 
simultaneously starting with expansions (88) to (90) in which use is made 
of the correct values of b^ and c Q obtained in the previous section. 

More details about this method of linear superposition of specialized 
solutions can be found in Appendix A where it is applied to a spherical 
perturbation due to finite external pressure. 

The integration has been carried out for the first five harmonics, 
treating each as if it were alone i.e. the coefficient A^ (i = 1,2, 3, 4, 5) 
in series (45) is set equal to one and all the others zero. An actual 
power perturbation G (0) would include several coefficients different from 
zero and with numerical values distinct of one. The effect of such a 
perturbation can easily be obtained by a linear superposition of the 
results presented here. 

The results describing the wave front and the physical variables behind 
it are presented in Table I; each harmonic is treated separately. As a 
reminder, the five first Legendre polynomials are given in the second 
column of Table I. The magnitude of the shape perturbation is given by 

of H (0). All X^s are 
positive indicating that an increase of laser power creates an outward dis- 
placement of the wave front and hence an increase of its velocity D. 

For increasing values of k, the magnitude of Xj^ diminishes: this shows that 
the shape perturbation becomes less when the power is more evenly 


the coefficient Xj^ of the Legendre expansion 


linear combinations of the specialized solutions ^ 
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distributed. The values of radial velocity, tangential velocity, density 
and pressure behind the absorption front are readily obtained from 
equations (51) to (54) combined with expansions in Legendre polynomials 
(95), (97), and (98) to (101). They read: 

V d = 1 - f (0) + e p k (cos 0) [-KT^ + (i _ f(°> + 5K) 3^] (131a) 

' W d = - « (1 - f (0) ) \ (131b) 

R d = ■ e p k (cos e) yr § tz (\ - 5 ^) (132a) 

P = 1 - + e p k (cos9 ) [-KA^ +{ 2 (1 - f^ o) ) 

+ 5 K } \] (132b) 

dP k 

Table I gives the numerical values of the coefficients of P, or when 

k d9 

A^ is set equal to one for each individual harmonic. The density perturba- 
tion corresponding to a power increase is negative: this simply results 
from the fact that the influence of the shape perturbation H (0) is 
greater than that of the power perturbation G (0) in all cases. As a 
result of the enhanced front velocity, the velocities and pressure perturba- 
tions at the boundary appear first positive, then negative. It must be 
pointed out that the absolute values of the perturbations of radial 
velocity, density and pressure increase with k; in fact they tend to 
asymptotic values for k ** °° as will be shown in section 3-4. The 
tangential velocity decreases with k as it should for a flow tending towards 
spherical symmetry. 

The radial dependence of velocities, density and pressure are plotted 
on Figures 4,5,6 and 7 for harmonics 1 to 4. The 0 - dependence is 
obtained by multiplying by the respective Legendre polynomials. Note 
that the radial variable A used here is not the ratio of the current 
radius r to the wave front radius r d , and consequently the wave front does 
not correspond exactly to 1 =1. The choice of A was made so that the 
radial and tangential variables A and 0 remain independent of each other. 
The tangential velocity remains negative as was imposed by the boundary 
condition that the velocity vector be perpendicular to the wave front. Its 
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magnitude decreases over the whole range of A for increasing k. The 
density perturbation exhibits a change of sign from negative to positive 
around A x .3. 

Although the character of the singularity near the origin changes 
between harmonic three and four (the velocities become bounded from 
harmonic four and up, whereas they were unbounded for the first three) no 
drastic change in the profiles appears to take place. This is discussed 

in more detail in Chapter III, Section 1. 

Density maps are presented in Figures 8 and 9; for the first two 
harmonics the constant density lines are shown. Compared to the spherical 
spark the mass repartition has been altered as follows: near the focus, 
an increase of density corresponds to enhanced power addition. In the main 
body of the flow, on the contrary, a decrease of density is associated 
with a positive power perturbation. Separating these two regions are 
circles of radius A = 0.34 and A = 0.3 for respectively k = 1 and k = 2. 
Near the boundary the density perturbation again follows the power 
perturbation. In the asymmetric case (k = 1) a net mass transfer appears 

to take place across the focal region from 6 = it towards 0=0. It 

should be remembered, however, that the small perturbation scheme breaks 
down near the origin so that the validity of what happens there is 
questionable. 

Figures 10 and 11 show velocity maps for k = 1 and k = 2. The vectors 
represent direction and magnitude of the velocity at different points of 
the field. Lines tangent to the velocity vectors are streamlines, not 
particle paths since the flow is unsteady. Physically, the vectors are 
proportional to the traces which would be left by flow particles made 
visible by flow visualization techniques on a short exposure time photo- 
graph. The bending of streamlines is determined by the distorted shape of 
the spark since velocity has to remain perpendicular to the wave front. 

In the asymmetric case, k = 1, there appears to be a net flow across the 
focal region from 6 = 0 towards 0 = it ; the zero velocity point is 
displaced on the cylindrical symmetry axis towards higher laser powers. 
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This, of course, doe not happen in the case of symmetric perturbation 

k = 2. 


3.4. Limiting Case of Large Wave Numbers 

For increasingly large wave numbers k it is observed in table 1 that 
the numerical value of the first root m^ of the indicial equation is 
positive and increasing without bounds; the fourth root m^ is negative and 
its absolute value is also increasing boundlessly. On the contrary both 
the real part m^ and the imaginary part m^ of the complex conjugate roots 
tend to finite asymptotic values. These behaviors have been confirmed by 
calculating the roots corresponding to the next ten harmonics ( k = 10, 

11, ...20). The values of m^ and m^ k are displayed on Figure 8. 

The purpose of this section is to show that in the limit k °° the 
original fourth order system of ordinary differential equations (104) to 
(108) reduces down to a second order system. Furthermore, the behavior 
of the physical quantities near the focus A = 0 is entirely described by 


the asymptotic values of the complex conjugate roots which will be 
calculated exactly. 

Physically, for large wave numbers, all quantities i.e. the power 
perturbation, the shape perturbation and the physical variables oscillate 
rapidly in the angular direction 9 thus giving rise to large 0 - deriva- 
tives. This is expressed by: 


d 

^■g - [P^ (cos9 )] ~ 0(k) when k 00 

Terms + cot 9 containing second derivatives of Legendre Polynom- 

® 2 

ials in equations (62) and (65) become 0 (k ) in equations (104) and (107) 
and the whole equation (64) containing first derivatives of P^ (cos0 ) 
becomes 0 (k) . A new ordering of terms is called for in order for the 
expansions (55) to (58) in powers of e to remain valid for large values of 

k. The perturbation quantities V (1) , W (1) , R (1) , P (1) have to be 0 <*) 
or less. 

It is noticed that one should expect a singular perturbation problem 
when k-*’ 00 due to the existence of two independent small parameters e and 
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The new ordering of terms corresponding to large k is suggested by 

equations (104) to (107). V k and R k remain of 0 (1) but W k is set of 

0 (Hr ) thus making P. of 0(— i) through equation (106). New tangential 
k Z k Z 

velocity and pressure are defined: 


W fc (X) - -4 W k (A) (133) 


Substituting (133) and (134) in equations (104) to (107) and neglecting 
in first approximation terms of 0 (— ^-) one obtains: 


. ,>>’ + 2 R 


v; + (r 

k 


(o) 

r > \ 


k-H R 

X 


(o) 


k X k 

(V (o) + 2 \ m o 


w u + (V (o) -A) r£ + 


(V (o) _ X) V' + (J2li + V (o) ) V. - - 

K n K 


(o)’ _ 

— \ - 0 


yP (o) V' + (P (o) + 2y 


R<°) 

(o) lf+1 p(®) — 

) V - Y V- W, = 0 
k k A k 


(135) 

(136) 

(137) 


Equation (106) contains terms which are all 0 (— ) ; it may be used as a 

_ k 

decoupled relation to calculate P k - The system (135) to (137) is only 
second order as can easily be seen by eliminating ^ and W k . 

Consideration of large wave numbers thus brings about the important 
simplification of reducing the order of the system by two. This is 
further exploited by specializing equations (135) to (137) to the neighbor 
hood of the origin. Only the first term in the expansions of the 
coefficients is retained so as to obtain: 


V a ' 2 V + h aX a ~\ 


k+1 


b,X 


«- 3 w k + (*„-!> A + 3 a 0 \ = 0 


U„ - »A K * 


+ a_) v. - c — — i 
° k 1 b 2 


-oH- 3 _ 

\ = 0 


V ' + 2 Vk - iS+1 W. = 0 


*k T A IT ~TT k 

In complete analogy to what was done for the fourth order system, 


(138) 

(139) 

(140) 


solutions to (138) and (140) are sought of the form: 


V k " A lk A 
% = A 2k A 


nioo 


(141) 

(142) 
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+ a- 3 


(143) 



An indicial equation is obtained for . It is quadratic as it should 

for a second order system and reads after some manipulations in which it 
k+1 

is remarked that — *• 1 when k-> <=°: 

k 

- (ao -1) m 2 -M. + [a-2)a 0 + n=l +a ] = 0 (144) 

^ n 00 ° a 0 -1 n ° 

The solution is the pair of complex conjugate roots: 

^ = 0.227 ±i 0.400 (145) 

These values are plotted on Figure 12. They show excellent agreement with 

the asymptotic values obtained from the fourth order indicial equation when 

k gets large (k = 20) . This also proves that the real part of m^ is less 

than one. The two purely real roots m^ and m^ completely disappear in 

im. 

this limit; this looks reasonable since the contribution of 4 tend to 

w 

be very small and can thus be discarded and on the other hand 1 gets 
very large, so violating the small perturbation assumption and has to be 
discarded too. 

This approach provides a rational way to reduce the order of the 
system and will be used in Chapter III when dealing with the full non- 
linear equations (32) to (35). 

It is easily recognized that boundary conditions (109) to (111) 
cannot be satisfied by equations (135) to (137). This comes as no sur- 
prise since the order of the system is only two and the number of 
independent boundary conditions is three. Specifically, in view of the 
definition (133) of equation (109) implies that is 0 (_^_ ) . A 

contradiction then appears between the definition (134) of and equation 

( 111 ). 


This difficulty can be resolved by introducing a narrow region near 
the front of the wave where derivatives with respect to X become large. In 
fact the singular perturbation aspect of the expansions suggests such a 
boundary layer approach. The proper ordering in that region is to keep 
V and of 0 (1) compatible with the boundary conditions and to 

make and the thickness of the boundary layer of 0 (— ) . One defines: 
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w = ± w 


k k . 

= 1 - lx 

k 


(146) 

(147) 


This ordering makes of 0 (1) so that the expansion in series of e 

is still valid. Furthermore (146) implies X. of 0 (A) through boundary 

k 

condition (109). Physically this corresponds to the perturbation of the 

wave front being much smaller, — , then the imposed power perturbation. 

k 

Thus boundary conditions (108), (110), (111) take only into account, in 
first approximation, the power perturbation A^. The perturbation in shape 
is reflected solely in the boundary condition (109) on the tangential 

velocity W, . 

k . i r . 

Equations (104) to (109) are transformed as follows. The coefficients 
dependent on the first order physical variables are expanded in Taylor 


series 


F (o) (X) 


= F (o) (l) 


( 1 ) * 


(148) 


where F(°) is any V(o), r(°), p(°). Use is made of the zeroth order 

boundary conditions (66) to (68). Terms of 0 (k) are retained in equations 

(104) , (105) , (107) compared to terms of 0 (1) and terms of 0 (1) are 

retained in equation (106) compared to terms of 0 (— ) . One obtains: 

k 


-JL A + Js±i 1 5 - f (o) d = o 

f(o) dV k f<°) k dA 

. A . dh _ 0 
+ 7T - 0 


(149) 


(150) 


+ P k - 0 


(151) 


Ta . f (o)) -^k + k+1 Y(1 . f <»)) ^ . £ <o) A . 0 (152) 

The last three equations are combined into a single equation for P^, 
k+1 

remarking that ~^r~ 1 when k -► 00 . 


2 \ ■ 0 


(153) 


where 


B 2 = y 


Y (l-f (o) ) - f (o) 


(154) 
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8 is a positive number and its numerical value, using the results of the 


spherically symmetric solution is: 8 = 3.390 or 6 = 1.840. 

Boundary conditions are given by (108) to (111) with = ix k : 

V k (l) = - 10^ (155) 


w k (D = - (1 - f (o) ) x k 


R k (l) = - 

8' 

p k a) = - 


(155) 

(156) 

(157) 


(158) 


The solution of (153) gives an increasing exponential and a decreasing 
exponential in A. The former is rejected since only bounded values of 
the physical variables are acceptable when A -*• °°, that is at the inner edge 
of the boundary layer. The solution of (149) to (151) and (153) subject 


to boundary conditions (155) to (158) is: 
_ 8 A 

V k = - e 


- ML e 
8 


Rk = - 


Yf(°) (l-f <o) ) 


- 1 ^ 
R 2 f<o) 2 


P k = - ^ e _ts A 

and i - + f^ 1 . 

k 6 l-f(o) 


(159) 

(160) 
(161) 


(162) 

(163) 


The validity of the boundary layer expansions near the front of the 
spark can be discussed in comparison with the numerical results obtained in 
Section 3.3. when k increases. Since calculations were carried out only 
for the first five harmonics, no accurate agreement should be expected but 
rather some general trends. 

First the boundary layer ordering implies that X k = 0 (^) ; thus for 
large values of k the expected behavior of subsequent X^s is given by 


. -i_ (i6 4) 

\ k+1 

In comparison this ratio is calculated for the first five harmonics using 


the values of X k obtained previously (equation (131). The result is shown 
in Table 2 and plotted on Figure 13. 
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k 

k 

\ + l 

k + 1 

calculated 

1 

0.5 

0.921 

2 

0.667 

0.904 

3 

0.750 

0.8996 

4 

0.800 

0.9015 


Table 2 


It is seen that the right trend appears in the numerical solutions when 
x k + 1 

k equals 3: the ratio — begins to increase and will presumably tend 

\ ; 

towards one. As an indication that k = 4 is still quite far away from the 

asymptotic limit k reference is made to Figure 12. 

Secondly, the boundary values obtained from (159) to (163) are: 

V 1 ) = P k (D = - 0.789; \(1) = - 2.71; W k (l)-> 0 (165) 

These are compared to those of the numerical solutions for the first five 
harmonics (see Figures 4 to 7). Again the correct trends emerge: v k (l) is 
negative with an increasing absolute value for increasing k's and the value 
(165) fork-*’ 00 fits in that pattern: see Figure 4; the same behavior is 
observed for R k (l) and P k (l). In contrast W k (l) is negative with a 
decreasing absolute value (Figure 5) which is also in agreement with (165). 

It may thus be concluded that the boundary layer expansions agree with 
the trends of the numerical solution for increasingly large k's. 

The boundary layer solution (159) to (162) must be compatible in the 
limit k-*°° with the solution of the inside equations (138) to (140). In 
particular and must become of 0 (A^-) as k -*■ 00 : this is certainly 

possible in view of the exponential decay (160), (162) of both and P^. 

It is noticed that R^ is finite at the inner edge of the boundary layer 
whereas is exponentially small; this of course does not imply that V k 
is small everywhere in the inner region. 

Matching of the boundary layer solution with the inner region 
solution thus appears possible. 
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CHAPTER III 


INVESTIGATION OF THE SINGULARITY AT THE FOCUS OF 
A NON-SPHERICAL SPARK 


1. Nature of the Singularity 

In Chapter II the behavior of all physical quantities was obtained 
near the focus X = 0. Equations (124) to (127) together with Table 1 
show that the two components of velocity and are oscillating and 
unbounded for the first three harmonics then become oscillating but of 
decreasing amplitude for subsequent harmonics. An analogous behavior is 
observed for the density except that the transition from unbounded 
amplitude to bounded amplitude occurs between harmonics two and three. 

The pressure P^ is oscillating and decreasing for all harmonics. 

Presented in this fashion the singularity at the origin may appear 
quite confusing. In order to get a better understanding of the physical 
situation, fluxes of basic properties such as mass, momentum and energy 
are computed across a sphere of small radius surrounding the focus. 

Since one of the basic assumptions made was that of non-viscous, non-heat 
conducting fluid, singularities of the type source or dipole are mathe- 
matically possible. In this problem, however, no source exists at the 
origin and such a singularity is not acceptable. This remark is used to 
provide a physical basis for the rejection of solutions depending on the 
fourth root of the indicial equation m^ . 

1.1. Integrated Flow Properties near the Focus 

Consider the fluxes of mass 7 ?? g > axial momentum'/^ and energy £ 

across a sphere defined by X = A with X << 1. This choice insures the 

o o 

reference sphere radius to be much smaller than the spark size at all 
times . 


vis 15 C °\ 
IB. ' C 

„ 2 


2 tt r^ sin 0 d0 = 2ir p n a\^ D ^ / R V X^ sin0 d9 (166) 

o o o o 

2 

cos 0- u. sin®) u 2 n r sin 0 d0= 2 if P 
0 r o o 

t 4n-2 jit £ (y cos 0 _ w sin0) VX sin0 d 0 (167) 

o o 
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(168) 


£ = 


u 2 + u 2 

p t ^) u 2tt r 2 sin 8 d 8 = 2 irp n^A^ t^ n 2 
o *■ Y-l Pro o 


2 2 

,ti „ .V + W 1 P v 2 

/ R ( — r + — r -) VA 

o 2 Y-l R o 


sin 8 d 9 


The mass and momentum fluxes increase with increasing time (respectively 
as t ^ and t ^ ) whereas the energy flux remains constant. This, of 
course, is a consequence of the constant power assumption. 

The 8 - integrals in equations (166) to (168) give the dependence 
on A^. In the mass flux m s the product RV contains successive powers of 
e : 


RV = R (o) V (o) + e (R (o) V (1) + V ( 0 ) R (1) ) + e 2 (R ( °M 2) + V (o) R (2) +' 


R (i) v ( 1 ) j + 0 (e 3) 


(169) 


Second order perturbations (superscript (2) ) are treated in Appendix B; 
their behavior is described in the neighborhood of the origin by equations 
(B 13) to (B 16). 


The integral from 0 to tt of the zeroth order term of (169) yields a 

mass flux whose dependence on A^ is: A^ ° + \ This simply corresponds to 

the spherically symmetric geometry. The 9 - integral of the 0 (e) term gives 

no contribution since the series of Legendre polynomials (98) to (101) 

start with £ = 1 and Legendre polynomials satisfy the integral relation (96). 
, 2 

There is, however, a contribution to 0 (e ) since the series (B 5) 
to (B 8) and the product R^^V^ contain the constant term P (cos 9 ) = 1. 


o 

If terms containing D's were retained in the first order expansions 
near the focus (124) to (127) , the leading terms of R^V^ as well as 
those of R^°^V^ ^ and V^°^R^ would be: A 2iDn + , thus contributing 


a mass flux: 


n. - k>< 


+ a-l 


(170) 


It is readily remarked that this always corresponds to a source of mass 
at the focus; even for the smallest possible value of m^, that is k = 1, 
the exponent of (170) is equal to - 5.091. 

Similarly it can be shown that the worst singularity m^ of the first 
order expansions also yields a source of energy at the focus. It does not 
give a source of linear momentum for the first harmonic although it does 
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for all higher harmonics. Since the physical problem does not allow for 
any kind of source, these considerations provide sufficient ground for the 
rejection of solutions based on D's in Chapter II, Section 3. 

It remains to show that the singularity related to the complex conju- 
gate roots does not create sources of mass, linear momentum or energy. 

2 

The leading contribution to 0 (e ) in (169) corresponds to the 

complex conjugate roots of harmonic 1: m ± im. and gives a mass flux 

r l L 1 

whose amplitude behaves as: 


MTtQ - A r 

S 0 1 


2m + a-1 


(171) 


The numerical value of the exponent is 0.840 which is physically acceptable. 

It may be remarked that if the leading harmonic happens to be higher than 

one, the exponent of X will always be greater than the above numerical 

o 1 

value, thus being a fortiori acceptable. 

2 

The linear momentum flux contains RV and RVW which expand respective- 

2 = R (o) (o) 2 + e (V (0) R (1) + 2R (0) V (0) V (1) ) + 0 ( e 2 ) (172) 

RVW = e r ( o) v (o) W (1) (173) 

The integral over 0 in (167) has no zeroth order part as it should for 
a spherically symmetric situation and the 0(e) contribution is, taking 


into account only the first harmonic 

2 „(o) 2 - , 4 „(o)„(o) /?7 


T/L * * C t 7 V 


R„ + T R 


Um ' ^ "q 1 3 ' "1 ' 3 “ ' N 'l "1' 


(V, + W n ) ] 


(174) 


This gives an oscillating momemtum flux the amplitude of which behaves as: 

AWi ) - eX V + a+1 (175) 

‘m o 1 

The numerical value of the exponent is 3.466, which is acceptable. 

The energy flux contains two parts: the flux of kinetic energy and 

the pressure work. Considerations analogous to what was done for 

arid }(l allow to show that there exists a zeroth order contribution and a 
m 

second order contribution. The latter has an amplitude behaving like: 

» 2 2m + a +1 . , 

(£) - e X q ^ a (176) 


The exponent is equal to 2.840. 

In summary the leading perturbations of the mass, momentum and energy 

2 2 

fluxes were shown to be respectively of 0 ( e ) , 0 (p ) and 0 (e ). 
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Furthermore, the singularity m^ gives rise to at least a source of mass 
and energy at the focus whereas the complex conjugate roots m + im. 
yield bounded fluxes to leading order. This justifies 
setting D's = 0 in expansions (124) to (127). 

1.2. Mathematical Description - Singular Perturbation 

The treatment will be limited to the mass flux , the other fluxes 

s 

being amenable to the same approach. Furthermore, only the leading 

harmonic of the first order perturbation, which will here be assumed to 

be harmonic one, is taken into account. Higher order perturbations 
3 4 

( e » e , *’) could be computed; it is clear that higher order equations 
will look like equations (B 1) to (B 4) of Appendix B with right-hand sides 
containing more complex functions of the lower order perturbations. The 
A - dependence of the amplitude of the third and fourth order quantities 
can be inferred to be: 


A <v <3> ) - x ^ 

.2(v <4) ) 

4m -3 
~ A r i 

(177) 

A <w <3) ) . x 3m tl " 2 


1 

. 4m -3 
-X t x 

(178) 

4(r < 3) ) - x 3m Cl + °- 5 

^t(R (4) ) 

4m + a-6 
- X "l 

(179) 

4<p <3) > - X 3 \ 1 + “' 3 

>C(P (4) ) 

4m + a-4 

> A r 1 , 

(180) 


The mass flux is then obtained from (166). The integral over 6 gives 

2 3 , 

terms of 0 (1), 0 (e ) , 0 ( e ) * • * The expansion has the form: 


rris K o 


, . «+l 2 2(m -1) f , , x . 3 3(m t1) .. 

(a)tca [ 1 + e A r f 0 (A)+eA r, 

ooo o!2o ol 


f, ( A ) + e 4 A 4<m r _1) f ( A ) +• • • * ] 
Jo o 1 4 o 


(181) 


where f , f and f are purely oscillating functions of A whose 
2 J 4 o 

arguments are respectively 2 i. In A, 3 in A and 4 m. In A. 

Similar series of powers of e can be obtained for the linear 
momentum and energy fluxes. 

Series (181) is uniformly convergent in Poincarre's sense as- long as 


. m -1 

t A r, <1 

q 1 

Thus there exists a region in the neighborhood of the origin of size 


(182) 
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1 


(183) 


0.615 

0(em-l)=0(e ) 

r l 

In which expansions (55) to (58) in powers of e break down. It must be 

remarked that the region of non-uniformity exists however small e is. 

The problem is thus a singluar perturbation problem (reference [52], 

Chapter V). It is interesting to observe that the usual warning sign of 

singular perturbation (reference [52]), namely e being the ratio of two 

independent lengths is not present in this particular case. 

Equation (182) indicates that for any m <1, there is a region of 

r k 

non-uniformity near the focus. Since it was shown (equation (145)) that 

m is always less than one the singular perturbation exists for all 
r k 

harmonics, including the limit k -+°°. 

2. Model Equations Near the Origin 

2.1. Attempt to use a Poincare-Lighthill-Kuo Method 

The usual Poincare-Lighthill-Kuo (P.L.K.) method, as presented for 
example in reference [52], Chapter 6, is not applicable to a system of 

four equations for four unknowns. This, of course, is due to the fact 

\ 

that a single straining of the independent variable cannot take care of 
the singularity of four dependent variables. What is needed is an individ- 
ual straining for each dependent variable. The proper strainings can be 
obtained using a method derived from Pritulo's remarks [53]. It requires 
the knowledge of the second order solution, as is done in Appendix B. 

The method is founded on the following consideration. A function 
F ( X,6 ) is expanded in series of e: 

F ( X, 0 ) = F (1) ( X, 0 ) + e F (2) ( X,0 ) + 0 (e 3 ) (184) 

F^ is singular near X = 0 and F^ even more singular. The dependent 
variables are also expanded in series of e : 

X =■ s + e X (s, ) + 0 (e 2 ) (185) 

0=<f>+e4> 1 (s, <J>)+0 (e 2 ) (186) 

Substituting (185) and (186) into (184) one collects powers of e after 
making use of the Taylor expansion of F^. The determination of the 
straining functions X^ and is such that the 0 (e) term of the expansion 
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in the new independent variables s, <p is no more singular than the 0 (1) 


term. This is written as: 


3F 


( 1 ) 


3X 


(s >* ) A i + a| 


(i) 


<s,4. ) .§ L + F (2) (s,<f> ) = '$ r 1 


where P* is any function of s, <j> no more singular than F 


( 1 ) 


(187) 

usually it 


is set equal to zero. 

In order to apply this method to the problem at hand, four functions 

t 

F are defined from the physical variables, by subtracting the regular 


zeroth order term: 

v * = — = v d) + e v (2) and w * = ” = w d) 

e e 

R* and P* are formed similarly. Using the results (124) to (127) of 


(188) 


Chapter II and (B 13) to (B 16) of Appendix B, individual relations (184) 
can be written for V*, W*, R* and P* from which individual strainings are 
obtained. 

However, no functions X^ and ^ could be found such that the mapping 
s,<P to X, 0 is one to one thus insuring functions V* to P* to be single- 
valued in X,0 . As an example consider V* and choose $ ^ to be zero; X^ 
is then obtained: 


X < Si 0 ) = - 


V (2) (s,6) 

( 1 ) 

3V 

lx (s ’ 0) 


(189) 


,(D 


Owing to the oscillating character of V and hence of its derivative ■, ' 

in both X and 0, X^ takes up infinite values of either sign, which in 
turn gives a non-acceptable transformation s to X. : 

Notice that this difficulty with a P.L.K. method is not specific to 
this particular problem but rather appears in any situation where the first 
order function F^ is oscillating. 

The method of matched asymptotic expansions could be used. 

As indicated by equation (183) an inner independent variable X is defined: 

= 0.615 (190) 


X =e X. with 0C ® ^ 


m -1 
r 

Expansions (124) to (127) show the ordering of the inner dependent variables 
to be: 

olX, 


V -c V; W »> ^ R = P = p<°> + e^^ 1 ) 


(191) 
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<p turns out to be constant and $*, are solution of a 

fourth order non-linear system of partial differential equations. This 
system was found to be no more tractable than the original non-linear 
system (32) to (35) so that this approach was interrupted. 


t 


2.2. First Order Non-linear Models 

In view of the difficulties encountered with the above standard 
techniques non-linear model equations are considered with the hope to 
determine if the singularity at the focus originates in the mathematical 
treatment of the problem or in the basic physical assumptions. 

The full non-linear equations (32) to (35) are used in the neighbor- 
hood of A = 0. For convenience they are slightly modified by extracting 
the leading part of the spherically symmetric solution near the origin: 

V = a A + v (192) 


W = w 
1 


= S = 


( £ A + S) 
x a-l b x 


(193) 

(194) 


P - c +A a-1 (c. A + p*) 

O 1 


(195) 

The system of non-linear partial differential equations for v, w, s, p* is 
derived in Appendix C, equations (C 6) to (C 10). It is used as the source 
of model equations. 

Notice that for small A ' s developing the physical quantities in powers 

of e is equivalent to linearizing system (C 6) to (C 10) . Model equations 

must thus be chosen so as to include non-linear terms. 

A first model is obtained by retaining only the linear and non-linear 

radial velocity terms in the r-momentum equation (C 8). After dividing through 

by (a + — -) and absorbing the multiplicative constant in a new definition 
J o n 

of v, one gets: 


with 


N Av + v = v v 

A A 


a -1 

N - o > 0 

4 - n ~ 1 

o n 


The numerical value of N is 11/6. 


The general solution of (196) is 


(196) 

(197) 
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obtained by inverting independent and dependent variables; with K an 
arbitrary constant it reads: 

A=Kv- N + 1 


N+l 


(198) 


If this solution is to be matched with the numerical solution of Chapter II 
developed in powers of £ outside the region of non-uniformity, v must 
be of 0 (£), Thus 

v = e v* (199) 

The solution for v*, now directly comparable to or is: 

*-N £ 

A = Kv + — v* (200) 

A schematic plot of v* (1) is presented in Figure 14 for positive K. 

The first term in (200) corresponds to the dotted line and indicates that 

v* goes to infinity when A goes to zero. This is the sort of behavior 

observed for the amplitude of in Chapter II. This, however, is strongly 

modified by the existence of the second term in (200) which corresponds 

to the non-linear right-hand side of (198) . It forces the solution to 

turn around before reaching A = 0; v* then tends to infinity along an 

N+l 

oblique asymptote: v* = — — A. Thus, however small £ ; there exists a 
forbidden region near A = o which the solution cannot enter; its size is 
given by: 

N 1 N_ 

A = e 1 ** 1 [K (N+1)] N+1 N N+1 (201) 

c 

N 

Notice that the characteristic power of e: = 0.647 is fairly close to 

the power X * 0.615 corresponding to the region of non-uniformity of the 
singular perturbation (see equation (190)). 

In this case it can be shown that the P.L.K. method gives exactly the 
answer (200) ; this is not surprising since the solution appears in the form 
of X being a function of v* expanded in powers of e. 

The non-linear term of (196), even though it is small, radically 
changes the character of the solution near A = o. Unfortunately, this 
change depends on the exact nature of the non-linear term. Suppose that 
instead of equation (196) another first order model equation is chosen: 


50 



( 202 ) 


N Av 


. 1 2 

+ v = r v 


This is quite plausible in view of the existence of a term 


w 


in 


equation (C 8) which very likely behaves as v /A . The solution of 
(202) is reduced to a quadrature at the transformation v = A z and 
£ = In A ; expressed in terms of v* it comes out to be, with K an arbitrary 
constant: 

A 

' ! V* =■ K (N+l) (203) 

N 

i 1 * 1 + K 


This is schematically plotted on Figure 14. The dotted line corresponds 

to equation (203) with £ set equal to zero; it is exactly the same as the 

dotted line corresponding to the first term of (200) . The non-linear term 

introduces a change near A = o which is quite different from the one 

observed with the first model equation: v* goes to zero with a slope of 
N+l 

— — thus exhibiting a typical boundary layer behavior. The boundary layer 
thickness , n defined as the value of A for which v^ is zero, comes out: 


N_ 

e N+l 


N 


A * e N+l (NK) N+l 
c 


(204) 


Notice that the same power of £ appears in both models. The P.L.K. method 
applied to equation (202) does not properly describe the boundary layer 
effect: instead of going to zero for ^ 0, v* goes to a finite value of 


The examination of these two non-linear first order model equations 
shows that the correct behavior of the solution near A = 0 is very sensitive 
to the exact nature of the non-linear terms retained. It is thus desirable 
to use a rational way to derive a model equation. 


2.3. Second Order Non-linear Model 

The second order model equation must have the feature that its linear 
part yields the oscillating behavior in A characteristic of V^. The 
pertinent non-linear terms can be rationally chosen in the limit of large 
wave number (k) . In that case, the linear system (104) to (107) reduces 
down from fourth order to second order. The non-linear system (C 6 ) to 
(C 10) is similarly reduced. 

For large k, the ordering of tangential velocity and pressure 
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perturbation is obtained by carrying over the results of the linear 
analysis (equations (133) and (134) ) to the non-linear system: 

w - 0 (£) (205) 

P* - 0 (~| ) (206) 

k 

Equations (C 6) to (C 10) are specialized to the axis of symmetry 6 = 0, 

ir on which the 0 derivatives of v and s vanish. Further, only equations in 

the neighborhood of A = 0 are of interest. For small A, provided that v 

a+1 


and w are leading over A 
0 


, equation (C 10) simply reads: 

♦ = 0 (207) 

The large k center line equations near A = o obtained from (C 7) and (C 8) 
are : 


1 VS 

- (o - 2) — v + (a -1) (As -s) - vs - (a -1) — = 0 
b^ o A A A 

(a -1) Av, + (a + — - ) v + a c, s + vv, = 0 
o Aon 1 A 

Elimination of s provides the second order equation in v: 

(a -2) aC l v + (a -l) 2 A 2 v + (a + — ) (a -l)(Av,-v) 
— — o AA o n o A 

1 2 2 2 
+ 2 (a -1) A vv + v v + (a -1) A v + v v , 
o AA AA o A A 

. . ,2 n 1 .12 . . . - n— 1. 

+ (3 a + ) v v - (a -1) v v - (a -1) (a + ) 

on A A on 


(208) 

(209) 


A ° (210) 


This equation contains a rather long series of non-linear terms all of 
which must be retained. It is seen from the ordering (190) , (191) of 
dependent and independent variables in the region of non-uniformity of the 
series of powers of e that all terms, linear and non-linear, in equation 
(210) are of 0 (e v. 

The explicit A- dependence in (210) can be eliminated by the 
transformation 

v = A z (211) 

€ - In A (212) 

With ' 3 the equation reads: 

2 d 

N z" + N z* + N z 

^^ 2^22 2 2 
- 2N^ z 2 " + z z" - z ' + zz' - z z' + z z ’ + M^z - 

- M, z 3 = 0 (213) 

4 
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The quantities N’s and M's are defined positive: 


N = (a -1) = .7333 
1 o 


M = 7 a - — — ~ = 1.800 
1 on 


N„ 


(a -1) (2 a = .8311 ; M = 4- a = .9091 

o o n 2 


N, => -(a-2) a (a - - ) = .4073 ; M = 4 a + — -( 3-1) 
3 o o n 3 on 

(a + — ) = 1.236 
o n 


(214) 


M = a-2 

4 


1.091 


In form (213) the second order equation for z is amenable to phase-plane 
techniques. See for instance reference [54]. With the standard notation: 

(215) 

(216) 


X 1 5 2 


X 2 = z' 
equation (213) becomes: 

dx 2 “ (X rV x 2 " M 2 (X 1 ' V (x i " V X 2 + M 4 x i (X l" R l )(X l +a o ) 


dX. 


dC 

dX. 


X 2 (X 1 “ N 1 ) 


1_ 

X„ 


(217) 

(218) 


The quantities R and R are defined as: 
N 1 

R = _3 = 1.3965 

a M,' 
o 4 


(219) 


= 1.2445 


N 1 M 2 


( 220 ) 


In principle, equation (217) is to be integrated in the phase plane 
X^, X2 taking singular points into account. The physical location * is then 
found by integrating equation (218) . 

For example, consider only the linear terms in equation (210); equation 
(217) simply reads: 


dX, 


N 


2 ^ - Vl 


n 2 x 

1 X 2 


( 221 ) 


It has one singular point (0,0) which is a focus, the characteristic roots 
of which are: - 0.773 ± 0.400i. The corresponding solution for v is: 


, .227 + i0.400 
v = X 


.227 - i0.400 


( 222 ) 


+ a 2 X 
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which is- in agreement with the results of Chapter II, Section 3.4. • 

The complete equation (217) exhibits five singular points: (0,0) 

which is a focus; S (R ,0), S (-a ,0) and S ( °°, ") which are saddle 
2 1 3 o 5 , 

points; and the singular point (N . ,“■)•• 

The local solutions in the neighborhood of S are given by : 




(223) 


The corresponding solutions for v are obtained from integration of equation 


(218): 


(v - X) 2 = X 2 (c 1 In X + c 2 ) (224) 

c^ and c 2 are arbitrary constants. 

This result is interpreted in the following way: singularity is 

associated with points v , X on the straight line v = N, X of the v, 

c c 1 

X plane. Consider the integral curves going through such a point; their 
behavior in the immediate neighborhood of v , X is obtained from equation 


(224) with: 


X = X + X 1 
c 

v = v + v ' 
c 

Expanding the In X for small X', one obtains: 


(225) 

(226) 


v' = + N X* + ✓cTTT r (227) 

1 1 c 

Thus the integral curves cross the line v = N^X with an infinite slope. 
Furthermore, there exists, on the axis of symmetry, a forbidden region near 

the origin: an integral curve coming from the half plane X > X turns back 

■ - c ' • 

toward that region after crossing v = NX . 

1 * *' ' . • .1 - ' . . 

The integral curves defined by equation (227) depend on two parameters: 
the position X^ at which the crossing occurs and the curvature (related to 
c^) of the square root function. 

The size of the forbidden region could, in principle, be determined 
by matching the small perturbation solution (222), depending on the constants 
and A^ to be determined from the boundary conditions, to the local non- 
linear solution (227) also depending on two free constants. As a first 
approximation of the relation of the size of the forbidden region toe , 
it may be assumed that the linear solution (222) is continued until it ! 
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) which is 


crosses the line v = N X. This yields a size X ^ 0 (e 

1 c 

in agreement with the ordering used in the singular perturbation approach: 

see equation (190) with m = 0.227. 

r 

In summary it is seen that a forbidden domain of approximate size 
1.291 

0 ( e ) occurs on the axis of symmetry. It seems likely that this 
conclusion also applies to other 9 — directions. 

3. Interpretation and Discussion 

The above analysis brings about the true nature of the singularity 
at the focus: rather than possessing unbounded oscillatory physical 
variables, as indicated by the linear theory, the flow displays a for- 
bidden region or hole near the origin. This, of course, is not physically 
acceptable. Since the mathematical treatment can no longer be criticized 
from the standpoint of singular perturbation, the source of the difficulty 
must be sought in the basic assumptions made to derive the flow equations. 
Self-similarity in this particular problem, besides other conditions, 

requires neglecting the counterpressure p as well as neglecting viscosity 

o 

and heat-conduction. The latter two are responsible for the infinite 

temperature and zero density obtained in the spherically symmetric geometry. 

It is shown in Appendix A that when counterpressure p is included in the 

o 

problem as a small perturbation t away from self-similarity, these features 
still prevail (see equations (A 37) and (A 39)). 

It is believed that the infinite temperature is the physical origin of 
difficulties appearing in non-spherically symmetric geometry. 

Each gas particle retains the entropy it got when it was processed by 
the narrow absorption layer. In particular, there exists always a point 
of zero velocity on the cylindrical symmetry axis, for example the focus 
itself if the power perturbation has a second symmetry axis at 0 = 90°. 
Consider the gas particles immediately adjacent to that point; they have a 
high entropy because they were processed at early time and have a different 
entropy according to the angle 9 at which they crossed the spark front. 

Very strong 9 - gradient of entropy (or temperature) are therefore present 
near the focus. Physically it seems that these must be smoothed out by 
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diffusion processes including radiative heat transfer. 

It is thus concluded that a proper physical approach of the flow 
near the focus should probably include momentum transfer and heat transfer 
effects. , 
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CHAPTER IV 


'• SUMMARY AND CONCLUSIONS 

Most experiments on laser breakdown of gases until now have relied 
upon a single focused beam, thus creating very asymmetric sparks and 
theories have been limited to one dimensional models. In this study, 
angular variations of power are introduced. Furthermore, from the point 
of view of energy concentration in a single spark it would be highly 
desirable to use a spherical array of lasers surrounding the plasma. 

This, of course, will provide a nearly spherically symmetric distribution 
of power to which the small perturbation analysis presented here is 
exactly applicable. 

Boundary conditions are derived for a narrow absorbing layer pre- 
ceded by a shock wave. The inner plasma is described as a non-viscous, 
non-heat conducting, and perfect gas. The thermal radiation emitted by 
the hot plasma is neglected compared to the laser beam power input, thus 
making the inner flow isentropic. Self-similarity is possible provided 

the counterpressure can be neglected and assuming constant power addition 

3/5 

in time; the size of the spark then grows like t 

The angular dependence is introduced as a perturbation in the power 
addition, of magnitude e, a small quantity. The self-similar equations 
and boundary conditions are split into a zeroth order set and a first order 
set. The former describes a purely spherical spark. A numerical outward 
integration yields the velocity, density and pressure profiles which are 
comparable to those obtained by Champetier and Al. [40] and Wilson and 
Turcotte [41]. The temperature at the origin is infinite and the pressure 
bounded. The singularity at the origin is identified as depending on the 
self-similarity variable A according to ^34/11. 

The first order perturbation equations and boundary conditions are 
linear and can be decomposed in Fourier series of the angular variable 9 . 
For the first five harmonics an outward numerical integration yields the 
velocity components, density and pressure profiles in function of A. The 
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size of the perturbed spark is determined; it turns out to be much smaller 
than the power perturbation which corresponds to a redistribution of the 
added energy in the whole flow. The wave front is strengthened by 
enhanced power addition, thus making velocity, density and pressure 
perturbations positive behind the front. Maps are presented of the stream 
lines and constant density lines for the first two harmonics. In the case 
of harmonic one, a mass transfer occurs near the focus towards the 
increased power addition; similarly the stream lines are pulled in that 
direction. 

The case of large wave numbers is examined. The tangential velocity 
and pressure perturbations become small in most of the field and the 
system of differential equations describing the flow reduces down from 
fourth order to second order. A boundary layer forms near the wave front 
to accommodate the boundary conditions. 

The small perturbation solution is valid away from the focus of the 
laser in which neighborhood the physical variables exhibit a singular 
behavior: they oscillate with reducing wave length. Analytical expansions 
in the variable A show that the singularity does not, to leading order, 
yield a source of mass, momentum and energy. The expansion in powers of 
e is, however, a singular perturbation near the focus. 

Several first order and second order model equations are considered 
with the aim of including non-linear effects near the origin. In 
particular, for large wave numbers, a second order equation is rationally 
obtained which describes the velocity perturbation on the axis of symmetry 
of the spark. Solution of the first order to that equation suggests, 
and study of the second order equation confirms, the existence of a 
forbidden region near the focus. 

Further improvements of the analysis presented here would include a 
detailed derivation of the size of the forbidden region by properly matching 
the non-linear solution with the small perturbation linear solution, in 
the large wave-number limit. This, however, still does not completely 
describe the physical situation. The basic assumptions should be reviewed; 
self-similarity likely should be excluded to properly take into account 
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heat transfer during early phases of the spark history when the infinite 
temperature singularity is established. Experimental evidence indicates 
the existence, inside the spark, of a hot core in which radiative transfer 
effects are probably strong; this may suggest the use of a constant 
temperature region near the origin. 
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APPENDIX A 


SPHERICAL LASER SPARK WITH COUNTERPRESSURE 


General Equations and Boundary Conditions 


The inclusion of counterpressure p introduces a new parameter 

o 

containing mass, length and time; this invalidates the dimensional arguments 
leading to self-similarity (see Chapter I) . The governing flow equations 
are thus the full time dependent equations (19) to (22). The analysis is 


restricted to purely spherical geometry i.e. 


— = 0 and u ~ 0 and constant 

90 0 


power k - 0. Defining A as previously (eq. (25) and (26)) and a modified 
time variable t as 


T = Y 


/ . n-1.2 

(n A t ) 


and non dimensionalizing the physical quantities as 

u r = n A t n_1 V ( A.t) (A2) 

p = P q R ( A,t) (A3) 

p - p n 2 A 2 t 2 n 2 P ( X,t) (A4) 

o 

equations (19), (20) and (22) can be rewritten, after cancellation of 
explicit t-dependence 

(V -A) R + R V + 2 — - 2 t R =0 (A5) 

A A A n x 

(V -a) V + — V + ;P - 2 — x V =0 (A6) 

A n R A n t 

(V- A) P +2 — P + vPV + 2y — - 2 — x P - 0 (A7) 

An A A n t 

4/5 

Note that the new independent variable t (which goes as t ) can be 
interpreted as the inverse square of a fictitious Mach number based on the 
sound velocity in the undisturbed gas and the wave front velocity that 
would pertain to a self-similar spherical spark. It was found more 
convenient, for perturbation purposes, to use this fictitious Mach number 
rather than the Mach number based on the actual front velocity as was done 
by Sakurai [I, II] in his study of ordinary blast waves. 

The spark boundary is an unknown function of time that may be expressed 
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m (t) (A8) 

The wave velocity is obtained by derivation with respect to time: 


r = At 
d 


n-1 , 2n-2 

n A t ( m - — - — t m ) 

dm 


(A9) 


in which m' ■= — . Defining m, which depends only on the unknown function 
dT 


m as: 


_ , . 2n-2 , 

m (t) = m - t m 


(A10) 


the boundary conditions are obtained from (15) to (18) for constant power 

P n 

P and uniform heat flux J = __L : at ^ = At m ("O 

47T 


V d - m (t) ( 1-f ) 

R = 1/f 
d 

P d = ^ T + [m (T )] 2 (1-f) 


with 


f = 


1 , , I_ J_,2 Y +1 1 

Y+l [y ~2 y (1 2 } 3 2. 3 

m v m 2 "an m m 


(All) 

(A12) 

(A13) 

(A14) 


A. 2. Small perturbation Assumption 

A small perturbation scheme from the self-similar spherical solution 

(zeroth order) is proposed by considering t as the small parameter. This 

type of approach can be regarded as a late time solution for which the 

-2/5 

velocity of the self-similar wave (which goes as t ) has decreased 
sufficiently so that the undisturbed sound speed is no longer negligible 
with respect to it and the effect of non-infinite Mach number 1 /t has to 
be taken into account . 

The function m (t) is expanded in a MacLaurin series, the first term 
of which is only retained 

* d = m (t) = 1 + * 2 t 1 (A15) 

The density ratio f is expanded in powers of t making use of the fact that 
the radical in (A 14) is non-zero for the overdriven zeroth order wave. 


One arrives at: 


f = f (<3) + t (K + K 2 X 2 ) 


(A16) 


.(o) 


where f and K are defined as in Chapter II (equations (49) and (50)) and 
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is another constant solely dependent on the zeroth order wave speed as: 


k 2---v — h d-^r 1 ' 2 

4tt an 2iran 

The boundary conditions on velocity, density and pressure are then 


(A17) 


developed in powers of t . This naturally suggests expansions of the flow 


variables over the whole field of A: 

V (A, t) = V (o) (A) + t V (1) (A) 

R (A.t) = R (0) (A) + t R (1) (A) 

P (A.t) = P (0) (A) + t P (1> (A) 


(A18) 

(A19) 


(A20) 


Expansions (A 18) to (A 20) are substituted into the flow equations 
(A 5) to (A 7) and the boundary conditions (A 11) to (A 14) . Collecting 
the terms of 0 (1) yields the zeroth-order problem which has already been 
solved in Chapter II § 2. The 0 (t) terms give rise to a system of three 
ordinary differential equations with variable coefficients for , R^^ 

and P^. They read with* ■— ® 

R (o) V (1) ' + (R (o) ' + 2 ^ y ~) V (1) + (V (0) -A) R (1) ' + (V (0) ' + 2 . 

^ - 2 ) R (1) = 0 (A21) 

(V (0) - A) V (1) ' + (V (0) ' - — ) V (1) - P — .-u- R (1) + P (1) r = 0 (A22) 

n R ( ° )2 R (0) 

+(,«»* + 2Y V (1) + (V (o) -») P (1 > V(Y V <0) ' + 

, 2y ) I' <1> - 0 (A23) 

The 0 (t) terms are collected in the boundary conditions which are further 

transferred from the unknown location A = 1 + A, t to the known location 

d 2 

A = 1 by a Taylor expansion. One obtains: 

V (1) (1) = - K + [ (l-f (o) ) — - K - V (0) ’(l)] -A (A24) 

n 2 2 

E<1) U> [ ^T2 + E<0) ' (1> 1 *2 (A25) 

P (1) (1) = 7 - K + t(l-f (o) ) — 2 - K - P (o) ' (1) ] A (A26) 

1 n / 2 

A 2.1. Expansions near A= 0 


The system of differential equations (A 21) to (A 23) is to be integrated 
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subject to boundary conditions (A 24) to (A 26). By elimination of the un- 
known constant these reduce down to only two, rendering impractical 

an inward numerical integration from the boundary X = 1. Instead the 
three linearly independent solutions are sought in the neighborhood of 
the origin. One is to be discarded as corresponding to a physically in- 
acceptable singularity. A numerical integration is conducted outwards 
with the two remaining solutions so as to satisfy boundary conditions at 
X = 1. 


Keeping only the leading term in the expansions of the coefficients 
of equations (A 21) to (A 23) , the flow equations near the origin read: 


b X °' 2 V (1) ’ + b <xX a-3 V (1) + (a -1) X R (1) '+ (3 a - 2 — ) 
1 1 o on 

R (1) = 0 

(a -1). X V (1) ' + (a - S=i) V (1) - ^ X' a+3 R (1) + ±- X"“ +2 
o % on ,2 b 1 

b l 

(1) ’ 

p: 7 - o 

1 (l) 

(1) ' V 7 Q) ' (1) 

Yc V v J + 2 Yc + (a -1) X P V } + 3Y a P V 7 = 0 

O O A O O 

.CD 


(A27) 

(A28) 

(A29) 


Note that for P behaving like a power of X the last two terms of 

cc 

equation (A 29) are of 0 (X ) compared to the last term of (A 28). They 
may thus be dropped to leading order. Three linearly independent solutions 
are sought, each of the form: 


v (1) - A x A* 

R (1) =A 2 A* + °- 3 

V W - A A x+ »- 1 
3 


(A30) 

(A31) 

(A32) 


Substituting these into (A 27) to (A 29) gives, after cancellation of all X 
dependence, a system of three linear homogeneous algebraic equations for 
A^, A^, A^ . x is solution of the indicial equation obtained by setting the 
determinant of this system equal to zero: 
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(x + a) 

(a -1) x + a - ^ 

o on 


(a -1) (x-1) - 2 


a (a - -) 
o o n 


x + cu. 1 1 


0 (A33) 


x + 2 


The three roots are immediately found: 

x » - 2 ; x 2 = - ( o-l) ; x 3 = 1 + 2 


n-1 1 


(A34) 


The first root x = - 2 yields an infinite velocity at the origin as well 

2 

as an infinite density. It can be seen that to order t the origin is a 
point source of mass and energy; this feature is not acceptable and the 
solution corresponding to the root x^ = -2 is rejected by setting the 
A's equal to zero. The system of algebraic equations is solved for the A's 


corresponding to roots two and three. For x^, one obtains 
A^ = 0 ; A^ = 0 ; A^ = arbitrary 


(A35) 


and for x„, the solution is: 1 

3 a a _ — 

A 1 = 0 J A 2 = arbitrary ; A 3 = - x ° + ^ — A 2 (A36) 

Note that in both cases none of the physical quantities blow up at the 
origin. In order to increase the accuracy of the numerical solution, the 
first non-zero term of the expansions is sought. To that effect, the first 
two terms in the expansions of the coefficients of equations (A 21) to (A 23) 
are retained and it is easily seen that the next order terms of and 


_(1) . . , . x+a.x + 2 a-3 

R are respectively A , a 


Denote the coefficients of these 


terms by B's. Remarking that for y = 5/3 the numerical values of the 

23 31 

retained roots are = - — and x 3 = — the expansions read: 
v (1) - B 1 A + B 2 A 65,11 


(A37) 


- B 1 A 12/11 + A* X 32/11 

P (l> . A 1 + A 2 X 54/11 
3 3 


(A38) 


(A39) 


where superscript 1 refers to the root x^ and superscript 2 refers to the 
root x^. It is interesting to remark that the solution corresponding to the 
root x^ has exactly the same type of singularity as the zeroth order self- 
similar solution, namely linear velocity in A and constant pressure; it is 
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1 2 

noted however that the coefficients are different. Choosing and 

1 2 1 2 

as the arbitrary constants the four constants B^ , , A^ and A^ are 

solution of an homogeneous algebraic third order system; for y = 5/3 and 


using equations (88) to (90) one obtains: 

R 1 _ 111 1 

B 2 44 b l B 1 


and 


.1 15 „1 

A, - - —7 c B 
3 4 o 1 


2 = 2025 x 2175 

A 2 154 x 374 

2 2175 2 

A = c B 

3 374 o 1 


« »? 
o 1 


(A40) 

(A41) 

(A42) 

(A43) 


A 2.2. Numerical Integration and Results 

The system of six ordinary first order differential equations formed 
by the zeroth order system plus equations (A 21) to (A 23) is integrated 
starting from the origin with expansions (88) to (90) and (A 37) to (A 39). 
A fourth order Runge-Kutta scheme is used. The boundary conditions are 
obtained from equations (A 24) to (A 26) by elimination of the unknown 
constant A. 


2 
,(D 


(1) ~ X 1 


R (1) (1) - X 3 


P (1) (l) - X 5 


(A44) 


in which the X's are known constants: 

X = - K ; X = (l-f (0) ) — - K - V (o) '(1) 
l z n z 

X - - -S- • X - - R <0) ’ (1) 

3 c (o)2 * 4 -(o)2 R U) 


(A45) 

(A46) 


X c = - - K ; X, « (1-f (o) ) 2 — - K„ - P (o) ’ (1) 
5 y 6 n 2 


(A47) 


Use is made of the linear . character of the system (A 21) to (A. 23) to 

12 

determine the constants B and B . A first Integration is carried out with 

1 2 1 ( 1 ) ( 1 ) ( 1 ) 

B^ = 1 and B^ = 0, call that solution , P ; the numerical 

values of these quantities at the boundary A = 1 , of course, do not satisfy 

1 2 

the boundary conditions. A second integration is made for B^ = 0 and B^ - 1, 
call that solution v 2 ^> R 2^’ P 2^' The ful1 s °l ution which must 
satisfy the boundary conditions is expressed as a linear combination of these 
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(A48) 


two specialized solutions i.e.: 

V (1) = V™ 

R (1) - R (1 1 ) + R™ Bj (A49) 

p (D = p (l) g l + p (l) g 2 (aso) 

1 2 

The constants B^ and B^ are solution of the second order algebraic system 
obtained by applying boundary conditions (A 44) on the full solution 

,( 1 ) 


(A 48) to (A 50). It reads: 

,(1) 


V ( J } (1) R ( ^(l) V v ;'(l) R x ;'(l) , -- 

( — - — - — - — ) b +( — — ) b: = — - — 

x 2 x 4 1 x 2 x 4 1 x 2 x 4 


( 


v ( l)(l, p v<» (1 , 

xT 5—> B i +< -E— 


p( “a> j x 

— t \ B = — 

X 6 1 X 2 


(A51) 


(A52) 


As a check the system was integrated once more using the so obtained values 
1 2 

of B^ and B^ in expansions (A 37) to (A 39) and it was found that the 
boundary conditions were accurately satisfied. 

The velocity, density and pressure perturbation profiles are plotted 


on Figure A and A is obtained 
A 2 - - 0.1993 


(A53) 


A 2.3. Discussion and Conclusions 

The numerical value of A^ is negative indicating an inward displace- 
ment of the boundary of the spark as well as a weakening of the wave 
characterized by a slower front velocity. This fact is also reflected in 
the perturbations of the physical quantities just behind the front. It 
should be noted that this applies at the actual front A = 1 - 0.1993 T and 
not at the fictitious point A = 1 (which is not even inside the spark) 
where boundary conditions are applied in the analysis. The value of a 
perturbation physical variable at the front is easily obtained: 

F ( J } = F (1) (1) + F (0) '(l) • A 2 


(A54) 


The particle velocity is seen to be lowered behind the detonation front and 
so is the density as might be expected for a weakened wave. The pressure, 
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however, is higher showing that" the contribution of the added external 
static pressure overrides the lowering of the dynamic pressure (see for 
example eq. (17)). 

The perturbation velocity is negative over the whole range of X and 
its profile remains remarkably linear up to A = 0.5 which is not surprising 
in view of the expansion (A 37) the linear part of which is expected to 
dominate over the power law part with the high exponent 65/11. The 
perturbation density and pressure are positive and a pressure increase is 
present at the origin. Note that since the density goes to zero, the 
temperature remains infinite at the origin as it should because no diffusion 
mechanism has been provided in this perturbation scheme. 

Comparing these results with Sakurai's [I, II] for constant energy 
blast waves, it seems apparent that a considerable simplification in the 
analysis and numerical solution has been gained by the adoption of the 
perturbation parameter t, inverse of a fictitious Mach number, rather than 
the actual front Mach number. This is reflected in the fact that the 
unknown constant does not enter the differential equations but is only 
present in the boundary conditions, thereby avoiding 5 the difficulty of 
having to deal with unbounded functions in the integration procedure. 

A direct comparison of the physical quantities distribution with those 

of a blast wave presented by Sakurai would require expressing the dependent 

r 

variables in terms of a new independent variable, say x = — , where r is 

r d • d 

the actual front. Noting that x is related to X by x = X (1 - t X^) this 
could easily be done 

F (x) = F (0) (A) + t [F (1) (A) - F (0)? (X) X £ X ] (A55) 

Nevertheless, the trends can be qualitatively described without that trans- 
formation. The perturbation profiles for density and pressure are markedly 
less steep near the boundary than for a blast wave. This, of course, is 
connected with the remark made in Chapter II § 2 that a Newtonian layer is 
not possible when energy is deposited right behind the leading shock wave. 
The velocity profiles exhibit the same negative behavior throughout and 
the density profiles appear to have the same trend in both cases. The 
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difference seems more marked in., the pressure profiles: -at the center a 
decrease in pressure is obtained for blast waves rather than the increase 

r * 

present here and a diminution of the pressure gradient just behind the 

shock wave is reported by Sakurai; it is not so for constant power laser sparks. 

r, * t ' ■ 

The applicability of this perturbation solution including counter- 
pressure for laser sparks can be found in certain special cases. 

Parameter t becomes significant when the sound speed in the undisturbed 
gas is high, for example if the gas is preheated. This physical situation 
has been actually realized in an experiment conducted by Ahmad and Key in 
1969 [III]. Instead of focusing the laser beam in a cold gas they first 
created a primary laser detonation wave in helium at 8 atm, after a 2 P 
sec delay they focused the beam of a second laser on the hot gas just 
behind the front of the primary spark (which by then had degenerated into 
an ordinary blast wave) . A secondary spark was obtained which developed 
in the heated gas inside the blast wave. 

Another application is related to the possibility of maintaining a 
laser spark, once initial breakdown has taken place, by using much smaller 
powers than those commonly used in experiments and hence creating much 
slower moving fronts. Raizer discussed that situation for air in 1968 
[IV], the criterion being that the gas should remain sufficiently hot 
( -20,000° K) to be able to absorb the laser photons on a short distance. 

Front velocities of the order of a few km/sec seem possible, thus making 
the inverse Mach number non-negligible with respect to one. 
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APPENDIX B 


S OME REM ARKS ABOUT THE SECOND ORDER SOLUTION 

• S • . — - r 1 . - - ■ ■ — — ■ , 

B.l . Second Order Equations 

2 

Expansions (55) to (58) in powers of e are extended to 0 (e ). 

v 2 

Collecting terms in £ the equations for the second order perturbation 


quantities are: 

R<°V? + (R<°>' 


+ )V®4, W («' 2) + cot0 „ (2 >) + 

A , . . A (\ ^ 

/_\t tt'®/ / n \ /-in f-t \ 


(V (o> -A ) R< 2 > + (V<°>' + 2 f ) R <2 >= - [V (1) R (2 > + 

R (1) V? + J W (1) R ( 0 1} + | R^ X) (2 V (1) + W ( 0 1} + cot 0 W (1) )] (Bl) 


(V*°* -X) V (2) + ( v (o) + ) v (2) 

A n 


p(o)' 


(0)2 “ + 7^ ‘ X 

1 W U?2 _ 1_ R (l) p (l) + 
X r (o)2 K ^ X 


- [v^v^ + y 

A A 0 

P (Q) ' R d)2 , 

r(°)3 1 


(v <o>_ A) „«> + ( Vl' + ^ w <2» + i_ 1 ,«). . [v U V J> 


+ 1 w u) w a> + i v u) w a> . : 

A 0 A 

,Co)„(2) ^ ,,(2) R 


fe» ( t ^ 1 


+ cot 0 VP ') + 


(V {o) -X) P ( jp + (y v (0) ’ + 2y ^ 


+ 2 — ) P (2) 
n 


[V (1) P <2 > + Y P (1) V <2) ♦ i V <1) P ( 0 1) + i Y p (1> (2V (1) + 

A A A 0 A 

Wg + cot 0 W (1) ) ] (B4) 

The solution of this fourth order linear system is sought as the sum of a 
complementary solution plus a particular solution. Since the left-hand 
sides of this system are identical to equations (62) to (65) the complemen- 
tary solution can be treated in the same way as the first order solution 
(see Chapter II, Section 3). In particular its behavior near the origin 
is given by (124) to (127). 

B.2. Particular Solution in the Neighborhood of X = 0. 


Consider the right-hand sides of (Bl) to (B4) . Their 0 - dependence 
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involves products of Legendre polynomials. If N is the highest harmonic 

ttl 

present in the first order solution, developing the product of two N — 
order Legendre polynomials in series yields Legendre polynomials up to 
order 2 N. The particular solution is thus sought in the form: 


v (2) 

- 1 
$5° 

w< 2) 

„(2) 

I 

l=o 

2N 

R 

= E 
JL=0 

p( 2 > 

2N 
= E 
A=0 


P (cos 9 ) V (\) 
l *■ 

dl p £ (co«e)w A (*) 

P (cos 9 ) lj(x) 

P (cos6 ) P^X) 


(B5) 

(B6) 

(B7) 

(B8) 


In order to find the leading term in the expansions of V to P near 
X = 0, only the first term in the expansions of the coefficients is kept on 
the left-hand side and only the leading terms of the right-hand sides are 
retained. 

Two cases are considered. First suppose that the terms D's are kept 

in expansions (124) to (127). Table 1 indicates that the strongest 

singularity m/ corresponds to the highest harmonic. The leading term on 

2m 1 — 1 

the right-hand side of (B 2) is thus seen to be X N The particular 


solution then has a X - dependence: 


V - x 

\ 

2N 

_ 

2m’ -1 

W_ - X 

N 

2N 

*=. 

2m' -1 + o-3 

R 2i< ' A 

N 


2m' - 1 + a 

P - A 

N 

2N 


(B9) 

(BIO) 

(Bll) 

(B12) 


Secondly, suppose that D's are set equal to zero. The strongest first 

order singularity is now given by the first harmonic k = 1. Hence the 

X - dependence is given by (B9) to (B12) where in' is replaced by m + im 

N r 1 

It may be noticed that it is physically more acceptable to 
have the strongest singularity at the origin corresponding to the most 
asymmetric perturbation i.e. the lowest harmonic rather than the opposite. 
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The complete second order particular solution has been calculated for 

the case: D's = 0 and the first order solution limited to its first 

harmonic. It comes out with m = - 0.6253 and m, = 0.5500: 

r i 

V (2) = A 2 ™ 1 - 1 [cos 2 0 (X x X 2im i + X* x“ 2im i + + sin 2 9 


(xj x 2im i + x^* x ' 2im i + ) ] 

w (2) = x 2mr 1 sine cose (x 2 x 2im i + x 2 x ^ im i + £ 2 ) 


R (2) = x 2m r + “" 4 [cos 2 0 (X 3 A 2im i + X* A 
(X^ A 2im i + X’* A" 2im i+/*) ] 


-2im 

i + £ 2 ; 

-2im, . 2 % . .2 


i + /? 3 ) + sin^ 0 


P (2) = A 2m r + “~ 2 [cos 2 9 (X 4 A 2im i + X* A ' 2i ®i + j*,) + sin 2 0 

(X* A 2im i + X '* A- 2im i + /-) 

4 4 4 

The constants X’s and Z's are solution of a system of fourteen linear 


(B13) 

(B14) 

(B15) 

(B16) 


algebraic equations. The numerical values are 
.-3 


X 2 = (0.12+il.88)10 

X 2 * (1.05+10.17)10 

X 3 = (0.50-11.52)10 

X, = (0. 80+i0.04)10 
4 


-3 


-3 


-3 


Z 1 = 1.77 10~ 3 X^ = (-0.06 - 10.94) 10~ 3 Z^= -0.89 10~ 3 
Z 2 = 0.22 10' 3 X^ - (-0.91+i0.39) 10~ 3 Z^ - 2.09 10 _3 


Z„ = -2.34 10“ 3 X’ = (0.12-i0.56)10~ 3 
3 4 


Z’= -5.78 10 
4 


-3 


Z, =-4.07 10 
4 


-3 
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APPENDIX C 


SECOND FORM OF THE FLOW EQUATIONS 


The complete non-linear self-similar flow equations (32) to (35) 
are transformed by defining: 

1 " ; ' 

R " ? •' . . , 

V = a ^ + v 

W ■ w 

s -^r x + s> ' , 

. a-1 , * 

P = C Q (c^ ^ + p ) 

With the definition: 


^ + 2v + w Q + cot 9 w 


the system reads: 

- (a-2) “ v - ~ 
b l b l 


VS 

+ (a -1) (* s, - s) + vs. - ( ot-1) t— + 


(Cl) 

(C2) 

(C3) 

(C4) 

(C5) 

(C6) 


ws, 


0 


X X 
V X + 


s $ = 0 


1 * 

+ 7 3 p 8 ■ 0 




n-l. 


L 

+ 

) V + 0 c, s 

+ - 

O 


n 2 1 

1 



w 1 

* 

W 

v 9 

- — + I s [X 

P X 



n-l. 1 * 


L 

+ 

) w + 7 — p 

+ 

O 


n b, K 9 




1 



* 

3 

X 


1 

X 


1 * * 
(a -1) a c. X + y c — $ + a c,v + y c, <f + (a -1) [X p, + ( a-1) p ] 
o 1 o^a^l^l o X 

+ v p^ + (a-1) E-2L + w P 9 + y ~ P <f = 0 

In the neighborhood of the origin, equations (C 7) to (C 10) are 

linearized as follows: 

X v + 2v + w + cot 9w = 0 
X i 9 

- ( o-2) — v + (a -1) * s, - (a -1) s = 0 
b^ o a o 

(a -1) A v, + (a + ^~ ^) v+ctc. s + 7— [ * p , + ( a -1) p ] = 0 
o Aon 1 A 

( V 1) Xw X + + a o ) W + ^ P d = 0 


(C7) 

(C8) 

(C9) 

(CIO) 

(Cll) 

(C12) 

(C13) 

(C14) 
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(CX^ 
(.CX^ 

v ® 

* S itv9 (CX*» 

- - h- h 

w 2 <yn\ £ot 

-V-. 

.* - *» ” :;l « — - s ;: *» - * * tu * 

- 1 ' „ay*'“ tva 

,u >- 18 e * t,^'- 

„ ec i» llsi ® 8 

„i e »» 4 S5,e ° 
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Power distribution: g (9 ) = 7 ^- [ 1 + e P, (cos0) ] 

4ir k 

Wave front: h(9)=l+£P^ (cos 9 ) 

Radial velocity behind the wave: V = 1 - + e P, (cos9) C 

d k Vk 

Tangential velocity behind the wave: W. = -e ^kCcosd) C.„ 

d “ wk 

de 

Density behind the wave: R. = — ~-r +e P, (cos0 ) C , 

d f (o) k Rk 

Pressure behind the wave: P. = 1 - + e P, (cos0 ) C 

d k Pk 


k 

P, (x) 
k 

X k 

C Vk 

C Wk 

V 

°Pk 

1 


X 


.1903 

.0491 

.0875 

-.1297 

.1366 

2 

3_ 

x 2 -l 

2 


.1754 

-.0164 

.0806 

-•3325 

.0642 


5 

x 3 -3X 


.1589 

-.0891 

.0730 

-.557 

-.0161 

3 


2 



34 

1 

CO 

o 

X 2 + 3 

.1429 

-.1586 

.0656 

-.773 

-.0939 

4 


8 



63x 5 - 70 x 3 +15x 

.1288 

-.2217 

.0591 

-.965 

-.1625 

j 

8 


Note: X = cos 9 

f (o) - .54041 l-f (o) - .45959 - 1.8504 

f (o> 

Table I: Spark Shape and Values of Velocities, Density 

and Pressure Behind the Wave Front 
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Figure 2. Velocity, density and pressure profiles of a 
spherically symmetric spark when y = 5/3. 
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Figure 3. Velocity, density and pressure profiles of a 
spherically symmetric spark when T"='T. 
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Figure 4. Radial velocity perturbation profiles for harmonics 
1 to 4 . 
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Figure 9 . Density map for harmonic 2 when 
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Figure 11. Velocity map for harmonic 2 when 





1.0 


NUMERICAL SOL'N 




Figure 13. Comparison of numerical values of wave strength with 
predicted values, for large k. 
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